AFRL-RV-HA-TR-201 1-1006 


Frequency-Dependent  Characteristics  of  Regional 
Seismic  Phases:  Propagation  of  Pn  in  Western  China 


Robert  L.  Nowack 
Wang-Ping  Chen 


Purdue  University 
550  Stadium  Mall  Dr. 
West  Lafayette  IN  47907 


Final  Report 


31  Jan  2011 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  IS  UNLIMITED. 


AIR  FORCE  RESEARCH  LABORATORY 
Space  Vehicles  Directorate 
29  Randolph  Rd 

AIR  FORCE  MATERIEL  COMMAND 
HANSCOM  AFB,  MA  01731-3010 


NOTICES 


Using  Government  drawings,  specifications,  or  other  data  included  in  this  document  for 
any  purpose  other  than  Government  procurement  does  not  in  any  way  obligate  the  U.S. 
Government.  The  fact  that  the  Government  formulated  or  supplied  the  drawings, 
specifications,  or  other  data  does  not  license  the  holder  or  any  other  person  or  corporation; 
or  convey  any  rights  or  permission  to  manufacture,  use,  or  sell  any  patented  invention  that 
may  relate  to  them. 

This  report  was  cleared  for  public  release  and  is  available  to  the  general  public,  including 
foreign  nationals.  Qualified  requestors  may  obtain  copies  of  this  report  from  the  Defense 
Technical  Information  Center  (DTIC)  (http://www.dtic.mil).  All  others  should  apply  to  the 
National  Technical  Information  Service. 

AFRL-RV-HA-TR-20 1 1-1006  HAS  BEEN  REVIEWED  AND  IS  APPROVED  FOR 
PUBLICATION  IN  ACCORDANCE  WITH  ASSIGNED  DISTRIBUTION  STATEMENT. 


//signature// 


//signature// 


ROBERT  J.  RAISTRICK  JILL  PARKER,  Maj,  USAF,  Acting  Chief 

Contract  Manager  Battlespace  Surveillance  Innovation  Center 


This  report  is  published  in  the  interest  of  scientific  and  technical  information  exchange,  and 
its  publication  does  not  constitute  the  Government’s  approval  or  disapproval  of  its  ideas  or 
findings. 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and  maintaining  the 
data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing 
this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202- 
4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently 
valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 

1 .  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

31-01-2011  Final  Report 

3.  DATES  COVERED  (From  -  To) 

04-20-2008  to  12-31-2010 

4.  TITLE  AND  SUBTITLE 

Frequency-dependent  Characteristics  of  Regional  Seismic  Phases: 

Propagation  of  Pn  in  Western  China 

5a.  CONTRACT  NUMBER 

FA87 1 8-08-C-0025 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

62601F 

6.  AUTHOR(S) 

Robert  L.  Nowack1  and  Wang-Ping  Chen2 

5d.  PROJECT  NUMBER 

1010 

5e.  TASK  NUMBER 

SM 

5f.  WORK  UNIT  NUMBER 

A1 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Purdue  University 

550  Stadium  Mall  Dr. 

West  Lafayette,  IN  47907 

8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory 

29  Randolph  Rd. 

Hanscom  AFB,  MA  01731-3010 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 

AFRL/RVBYE 

11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 

AFRL-RV -H  A-TR-20 11-1006 

12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 


Approved  for  Public  Release;  Distribution  Unlimited. 


13.  SUPPLEMENTARY  NOTES 

*Dept  of  Earth  and  Atmospheric  Sciences,  Purdue  University,  West  Lafayette,  IN  47907 
2Dept  of  Geology,  University  of  Illinois-Urbana,  Urbana,  IL  61801 

14.  ABSTRACT 

The  propagation  of  Pn  waves  has  been  investigated  over  a  large  region  in  western  China  using  regional  events  recorded  by  the 
Hi-CLIMB  array.  Seismic  attributes,  including  arrival  times,  Hilbert  envelope  amplitudes,  and  pulse  frequencies,  are  modeled 
to  investigate  how  velocity  structures  in  the  crust  and  the  upper  mantle  affect  the  propagation  of  the  Pn  wave  in  Tibet.  Travel- 
times  are  first  modeled  using  3D  ray  tracing  to  construct  a  3D  velocity  model  beneath  the  Hi-CLIMB  array.  The  spectral 
element  method,  which  provides  complete  numerical  wave  solutions  in  laterally  varying  media,  was  then  used  to  model  the 
seismic  attribute  data.  Velocity  and  attenuation  models  beneath  the  Hi-CLIMB  array  are  constructed  and  trade-offs  between 
structural  effects  and  attenuation  are  investigated.  The  objective  of  the  study  was  to  obtain  self-consistent  models  of  Pn 
propagation  in  western  China  free  from  assumptions  such  as  a  frequency-independent  geometric  spreading.  To  this  end,  our 
results  should  advance  efforts  in  isolating  effects  of  frequency-dependent  propagation  from  those  of  pure-inelastic  attenuation 
(Q),  leading  to  improved  methodologies  for  discrimination  and  yield  estimates  at  regional  distances. 


15.  SUBJECT  TERMS 

Regional  Seismic  Wave  Propagation,  Pn  Waves,  Hi-CLIMB  array 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF  ABSTRACT 

18.  NUMBER 
OF  PAGES 

19a.  NAME  OF  RESPONSIBLE  PERSON 

Robert  Raistrick 

a.  REPORT 

UNC 

b.  ABSTRACT 

UNC 

c.  THIS  PAGE 

UNC 

SAR 

80 

19b.  TELEPHONE  NUMBER  (include  area 
code) 

781-377-3726 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  Z39.18 


Table  of  Contents 

1 .  Summary  1 

2.  Introduction  1 

3.  Travel-Time  Modeling  of  Regional  Events  in  Tibet  2 

4.  Spectral  Element  Modeling  of  the  Pn  Phase  6 

5.  Moment  Tensors  of  Regional  Events  in  Tibet  11 

6.  Specification  of  Moment  Tensors  for  2D  SEM  Modeling  14 

7.  Modeling  of  Pn  Seismic  Attributes  in  Tibet  15 

8.  Seismic  Profiling  with  Teleseismic  Waves  27 

9.  Conclusions  29 

10.  Recommendations  30 

References  3 1 

Appendix  A.  Calculation  of  Seismic  Attributes  35 

Appendix  B.  Validation  of  the  Spectral  Element  Method  for  Pn  Waves  37 

Appendix  C.  Preprint  of  Paper  by  Griffin  et  al.  (2011)  49 

List  of  Symbols,  Abbreviations,  and  Acronyms  72 


iii 


Figures 


1 .  Map  of  the  stations  of  the  Hi-CLIMB  array  in  Tibet  (triangles)  and  the  epicenters  of  the 

selected  regional  earthquake  sources  (circles).  Black  circles  are  the  hypocenters  reports  by  the 
PDE  catalogs  and  Engdahl  et  al.  (1998),  and  the  white  circles  are  the  epicenters  nearest  to  the 
array  that  have  been  relocated  in  this  study.  The  dashed  lines  indicate  the  locations  of  north-south 
slices  of  the  P-wave  velocities  from  travel-time  modeling,  with  the  East  Line  at  84  degrees,  the 
Central  Line  at  85  degrees,  and  the  West  Line  at  86  degrees  in  longitude.  2 

2.  Ray  diagrams  calculated  using  the  3D  complete  ray  tracing  code  (CRT)  for  the  Event  LIE 

located  by  the  dot.  3 

3.  These  plots  show  the  observed  first-arrival  travel  times  and  the  calculated  travel  times  (lines 

and  crosses)  for  the  relocated  nearby  events  using  the  average  ID  crustal  model  on  top  and  the 
refined  crustal  velocity  model  on  the  bottom.  4 

4.  These  plots  shows  observed  P-wav  e  travel  times  (circles)  and  calculated  travel  times  (1  ines 
and  crosses)  for  the  far-regional  events,  (A)  Lie,  (B)  Lla,  (C)  Lli,  (D)  Llj,  (E)  Llk,  and  (F)  E3b. 
The  horizontal  bars  highlight  the  locations  of  the  transition  for  shoaling  of  the  Moho. 

5 

5.  The  3D  Moho  and  upper-mantle  velocity  structure  found  in  this  study  for  the  region  of  the  Hi- 

CLIMB  array,  with  the  left  plot  showing  a  view  from  the  west  and  the  right  plot  showing  a  view 
from  the  south.  5 

6.  North-south  cross-sections  of  the  3D  velocity  model  of  Griffin  et  al.  (2011)  along  longitudes 

of  84,  85  and  86  degrees.  6 

7.  The  top  plot  shows  the  P-wave  vel  ocities  for  a  si  mple  layer  over  a  gradient  with  an  upper- 
mantle  velocity  gradient  of  .004  1/s.  The  bottom  plot  shows  a  simplified  P-wave  velocity  model 
from  travel-time  modeling  of  regional  events  from  Hi-CLIMB,  showing  a  step  in  the  Moho. 

7 

8.  SEM  synthetic  P-wave  A)  travel-times,  B)  amplitudes  and  C)  pulse  centroid  frequencies  for  a 

flat  Moho  with  several  velocity  models  with  different  upper-mantle  velocity  gradients  of  0.0, 
0.002,  0.004,  and  0.006  1/s.  For  this  case  no  Q  is  included.  8 

9.  SEM  synthetic  P-wave  A)  travel-times,  B)  amplitudes  and  C)  pulse  centroid  frequencies  for  a 

flat  Moho  with  several  velocity  models  with  different  upper-mantle  velocity  gradients  of  0.0, 
0.002,  0.004,  and  0.006  1/s.  For  this  case  a  Q  of  150  is  assumed  in  the  crust  and  400  in  the 
mantle.  9 

10.  SEM  synthetic  A)  travel-times,  B)  amplitudes  and  C)  pulse  centroid  frequencies  for  a  step  in 

the  Moho,  as  shown  for  t  he  right  plot  in  Figure  7.  The  crosses  for  source  B  near  1000  km  in 
Figure  7  (right  plot)  and  recorded  to  the  south,  and  the  circles  are  for  a  source  A  near  0  km  and 
recorded  to  the  north  in  Figure  7  (right  plot).  10 

1 1 .  SEM  synthetic  A)  travel-times,  B)  amplitudes  and  C)  pulse  centroid  frequencies  for  a  step  in 
the  Moho,  as  shown  for  the  right  plot  in  Figure  7.  The  crosses  for  source  B  near  1000  km  in 
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Figure  7  (right  plot)  and  recorded  to  the  south,  and  the  circles  are  for  a  source  A  near  0  km  and 
recorded  to  the  north  in  Figure  7  (right  plot).  1 1 

12.  Inversion  of  three-component  waveforms  for  the  focal-mechanism  of  event  E3b  using  data 
from  the  Hi-CLIMB  array.  The  top  right  plot  shows  the  data  misfit  with  a  narrow  minimum  at  6 
km  and  and  a  broad  flat  minimum  at  19  km,  and  the  upper  left  plot  shows  the  19  km  solution. 

The  lower  plots  show  the  matches  between  the  observed  and  calculated  waveform  data.  12 

13.  Focal  mechanism  solutions  for  Event  E3b.  a)  A  comparison  global  CMT  focal  mechanism 

for  an  earthquake  located  about  70  km  distant  from  Event  E3b.  b)  The  focal  mechanism  solution 
for  Event  E3b  with  a  data  misfit  minimum  at  a  depth  of  6  km.  c)  The  focal  mechanism  solution 
for  Event  E3b  with  a  dat  a  misfit  minimum  at  a  depth  of  19  km.  d)  The  focal  mechanism  for 
Event  E3b  for  a  shift  in  epicenter  of  10  km  to  the  east  resulting  in  a  least  misfit  minimum  at  the 
shallower  depth  of  7  km.  All  solutions  indicate  East-West  extension.  13 

14.  The  left  plot  shows  the  3D  focal  mechanism  for  a  particular  double-couple  moment  tensor 

source  with  four  angles  given  by  dashed  lines.  The  Pg  and  Pn  take-off  angles  are  shown  by  small 
circles.  The  right  plot  shows  the  amplitudes  on  the  focal  sphere  of  the  2D  m  oment  tensor  for  the 
azimuthal  angles  given  in  the  left  plot  as  a  function  of  take-off  angle  derived  from  the  3D 
moment  tensor.  1 4 

15.  The  left  plot  shows  a  focal  mechanism  with  four  azimuthal  angles  given  by  dashed  lines.  The 

Pn  and  Pg  take-off  angles  are  shown  by  the  small  circles.  The  right  plot  shows  the  SEM  envelope 
amplitudes  for  the  azim  uthal  angles  given  on  th  e  left.  Different  symbols  represent  different 
azimuthal  angles .  1 5 

16.  The  left  map  shows  the  original  and  slightly  revised  focal  mechanism  for  event  LIE.  The 

middle  map  shows  the  focal  mechanism  TL-1  and  an  alternative  focal  mechanism  TL-2  for  event 
E3B.  The  right  map  shows  the  focal  mechanisms  for  the  cluster  of  south  events  LI  A,  L1I,  LI  J, 
and  L1K.  The  triangles  show  the  Hi-Climb  stations,  where  the  solid  triangles  along  with  the  solid 
lines  show  the  orientation  angles  used  for  the  specification  of  the  moment  tensor.  The  dashed 
lines  show  the  separation  of  ranges  used  to  specify  the  moment  tensor.  1 5 

17.  This  map  shows  the  location  where  2D  P-wave  velocity  slices  were  obtained  from  3D 

velocity  model  of  Griffin  et  al.  (201 1).  Triangles  represent  the  Hi-Climb  array.  The  circles  show 
the  events  used  by  Griffin  et  al.  (2011)  for  the  3D  velocity  model  where  the  open  circles  are 
relocations.  The  dashed  lines  show  the  West,  Central,  and  East  lines  for  the  2D  slices.  16 

18.  The  plots  show  the  P-velocities  for  3  different  slices  of  the  3D  model  of  Griffin  et  al.  (201 1) 

which  constant  upper  mantle  velocity  gradients  but  different  upper  mantle  top  velocities.  The 
small  circles  show  location  of  the  events  used  with  each  velocity  slice.  17 

19.  These  plots  show  the  calculated  travel-times  of  the  selected  events  with  no  attenuation  in  the 

SEM  models  compared  with  the  observed  travel-times.  Circles  are  the  SEM  calculated  travel- 
times  and  crosses  are  observed  travel-times.  1 8 

20.  These  plots  show  the  peak  envelope  amplitudes  of  the  selected  events  with  no  attenuation  in 
the  SEM  models  compared  with  the  observed  envelope  amplitudes.  The  circles  are  the  SEM 
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calculated  values  and  the  crosses  are  observed  envelope  amplitudes.  The  distance  ranges  for  the 
specified  orientation  angles  for  the  moment  tensor  are  also  shown.  For  event  LIE,  the  modified 
moment  tensor  is  used  and  for  event  E3B  the  TL-1  moment  tensor  is  shown.  19 

21 .  These  plot  show  the  calculated  instantaneous  pulse  frequencies  of  the  selected  events  with  no 

attenuation  for  the  SEM  models  compared  with  the  observed  instantaneous  and  pulse  frequencies. 
The  circles  are  the  SEM  calculations  and  the  crosses  and  squares  are  the  observed  instantaneous 
and  centroid  frequencies.  For  event  LIE,  the  modified  moment  tensor  is  used  and  for  event  E3B 
the  TL-1  moment  tensor  is  shown.  20 

22.  These  plots  show  variable  upper  mantle  Q  models  for  each  slice  determined  from  modeling 
seismic  attributes  from  the  selected  regional  events  assuming  the  velocity  slices  in  Figure  18.21 

23.  These  plot  show  the  SEM  calculated  peak  envelope  amplitudes  of  the  selected  events  with 

the  attenuation  models  shown  in  Figure  22  and  the  velocity  slices  given  in  Figure  18  compared 
with  the  observed  envelope  amplitudes.  The  circles  are  the  SEM  calculated  values  and  the  crosses 
are  observed  envelope  amplitudes.  The  distance  ranges  for  the  specified  orientation  angles  for  the 
moment  tensor  are  also  shown.  For  event  LIE,  the  modified  moment  tensor  is  used  and  for  event 
E3B  the  TL-1  moment  tensor  is  shown.  22 

24.  These  plot  show  the  SEM  calculated  pulse  frequencies  of  the  selected  events  with  the 

attenuation  models  shown  in  Figure  22  and  the  velocity  slices  given  in  Figure  18  compared  with 
the  observed  pulse  frqeuencies.  The  circles  are  the  SEM  calculations  and  the  crosses  and  squares 
are  the  observed  instantaneous  and  centroid  frequencies.  For  event  LIE,  the  modified  moment 
tensor  is  used  and  for  event  E3B  the  TL-1  moment  tensor  is  shown.  23 

25.  The  plots  show  the  P-velocities  for  3  different  slices  of  the  3D  model  of  Griffin  et  al.  (201 1) 

which  variable  upper  mantle  velocity  gradients.  The  small  circles  show  location  of  the  events 
used  with  each  velocity  slice.  24 

26.  The  plots  shows  constant  upper  mantle  attenuation  used  with  the  velocity  models  in  Figure 

25  with  a  laterally  variable  upper  mantle  velocity  gradient  between  the  Lhasa  and  Qiangtang 
terranes.  25 

27.  These  plots  show  the  peak  envelope  amplitudes  of  the  selected  events  with  a  constant 

attenuation  in  the  upper  mantle  for  the  SEM  models  compared  with  the  observed  envelope 
amplitudes.  The  circles  are  the  SEM  calculations  and  the  crosses  are  observe  envelope 
amplitudes.  26 

28.  These  plots  show  the  instantaneous  pulse  frequencies  of  the  selected  events  with  a  constant 

attenuation  and  a  variable  upper  mantle  velocity  gradient  for  the  SEM  models  compared  with  the 
observed  instantaneous  and  centroid  frequencies.  The  circles  are  the  SEM  calculations  and  the 
crosses  and  squares  are  observed  instantaneous  and  centroid  frequencies.  27 

29.  North-south  trending  profile  across  southern  and  central  Tibet,  showing  topography  on  top,  a 

P-wave  reflection  profile  constructed  using  the  phase  SsPmp,  and  the  estimated  crustal  thickness 
(Panels  a,  b,  c,  respectively).  28 
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30.  A  comparison  between  crustal  thickness  estimated  from  wide-angle  P-wave  reflections 
(dots)  from  the  teleseismic  SsPmp  phase  and  an  image  of  the  Tibetan  lithosphere  obtained  from 
Gaussian  beam  migration  of  direct  P-  to  S-wave  conversions  (from  Tseng  et  al.,  2009;  and  the 
Gaussian  beam  migration  from  Nowack  et  al.,  2010).  29 

A. l.  This  shows  an  example  of  a  data  trace  from  the  Hi-CLIMB  array  with  a  window  function  at 

the  top  and  a  windowed  trace  and  envelope  at  the  bottom.  32 

B. l.  Source  time  functions  for  the  reflectivity  calculation  (solid  line)  with  t  =  0.2  sec  and  the 

SEM  calculation  with  a0=7A  (dashed  line).  34 

B.2.  SEM  synthetics  of  model  Moho-1  with  no  3D  correction  (top)  and  SEM  synthetics  of  model 
Moho-1  with  3D  distance  and  pulse  shape  corrections  (bottom).  Top  and  bottom  amplitudes  are 
scaled  by  the  maximum  amplitude  of  all  traces.  39 

B.3.  A  comparison  between  the  SEM  calculations  (dashed)  and  the  reflectivity  results  (solid)  of 
Braile  and  Smith  (1975)  for  a  constant  velocity  30  km  crust  over  a  constant  velocity  mantle.  The 
amplitudes  of  the  reflectivity  results  and  the  3D  corrected  SEM  results  are  both  multiplied  by 
distance  in  the  plot.  40 

B.4.  A  comparison  between  the  results  from  the  SEM  calculations  (dashed)  and  the  reflectivity 
results  (solid)  of  Braile  and  Smith  (1975)  for  a  constant  velocity  30  km  crust  over  a  mantle  with  a 
velocity  gradient  in  the  top  5  km.  Amplitudes  of  the  reflectivity  results  and  the  corrected  SEM 
results  are  both  multiplied  by  distance  in  the  plot.  41 

B.5.  This  shows  a  co  mparison  between  the  SEM  calculations  (crosses)  and  the  asy  mptotic 
calculations  (background  plot)  of  Cerveny  and  Ravindra  (1971)  for  the  pure  and  interference  head 
waves.  C+  wave  is  the  interference  head  wave,  C  0  is  the  diving  wave,  and  th  e  131  wave  is  the 
pure  head  wave  case  (Modified  background  plot  from  Cerveny  and  Ravindra,  1971).  42 

B.6.  This  shows  amplitude  spectra  as  a  function  of  distance  for  the  interference  head  wave  case. 
The  thin  black  curves  are  from  Cerveny  and  Ravindra  (1971)  and  the  thicker  curves  are  from 
SEM  calculations  for  the  same  model  (Background  plot  from  Cerveny  and  Ravindra,  1971).  43 

B.7.  This  shows  the  comparison  of  the  SEM  traces  with  the  analytical  calculations  for  a 
homogeneous  viscoelastic  model  with  a  Qp  of  1 12.5.  44 

B.8.  This  shows  the  envelope  amplitude  of  the  SEM  traces  with  the  envelope  amplitude  of  the 
analytically  calculated  traces  given  in  figure  B.7.  45 

B.9.  This  shows  a  comparison  of  the  frequency  spectra  of  the  SEM  traces  with  analytical 
calculations  for  the  traces  given  in  figure  B.7. 
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1.  SUMMARY 


The  propagation  of  Pn  waves  has  been  investigated  over  a  large  reg  ion  in  western 
China  using  regional  events  recorded  by  th  e  Hi-CLIMB  array.  S  eismic  attributes, 
including  arrival  times,  Hilbert  envelope  amplitudes,  and  pulse  frequencies,  are  modeled 
to  investigate  how  velocity  structures  in  the  crust  and  the  upper  m  antle  affect  the 
propagation  of  the  Pn  wave  in  Tibet.  Travel-times  are  first  modeled  using  3D  ray  tracing 
to  construct  a  3D  velocity  model  beneath  the  Hi-CLIMB  array.  The  spectral  elem  ent 
method,  which  provides  com  plete  numerical  wave  solutions  in  late  rally  varying  media, 
was  then  u  sed  to  m  odel  the  seismic  attribute  data.  Velocity  and  a  ttenuation  models 
beneath  the  Hi-CLIMB  array  are  constructed  and  trade-offs  between  structural  effects 
and  attenuation  are  investigated.  The  objective  of  the  study  was  to  obtain  self-consistent 
models  of  Pn  propagation  in  western  China  free  fr  om  assumptions  such  as  a  frequen  cy- 
independent  geometric  spreading.  To  this  end,  our  results  should  advance  efforts  in 
isolating  effects  of  ffeque  ncy-dependent  propagation  fr  om  those  of  pure-inelastic 
attenuation  (Q),  leading  to  im  proved  methodologies  for  discri  mination  and  yield 
estimates  at  regional  distances. 

2.  INTRODUCTION 

Using  seismic  data  from  regional  earthqua  kes  recorded  by  the  Hi-CLIMB  array, 
we  have  investigated  trade-offs  between  velo  city  structure  and  attenuation  in  Tibet.  Of 
more  than  30  high-quality,  regional  seism  ic  profiles  that  have  been  assembled,  we  have 
selected  10  which  show  excellent  crustal  and  P„  arrivals  for  further  an  alysis.  We  first 
model  P-wave  arrival-times  to  constrain  the  velocity  structure  in  the  crust  and  the  upper 
mantle  in  central  and  western  Tibet.  Travel-times  from  four  events  along  the  Hi-CLIMB 
array  provide  details  on  crustal  velocities,  and  six  events  at  regional  distances  to  the  array 
provide  further  constraints  on  Moho  structure  and  upper  mantle  lid  velocities.  We  use  3D 
ray  tracing  to  m  odel  the  travel-times,  and  the  results  indicate  that  both  Moho  structure 
and  mantle-lid  velocities  in  the  region  are  three-dimensional  in  nature  but  approximately 
follow  the  trend  of  the  Bangong-Nujiang  Sutu  re  (BNS).  Although  only  a  lim  ited  number 
of  events  were  used  for  the  travel-time  modeling,  the  results  are  consistent  with  earlier 
results  from  teleseismic  imaging  using  the  Hi-CLIMB  array. 

Attributes  of  the  Pn  waves,  including  envelope  amplitudes  and  instantaneous 
frequencies  (see  Appendix  A)  ,  are  next  m  odeled  using  th  e  spectral  elem  ent  method 
(SEM)  using  regional  earthquakes  recorded  by  the  Hi-Climb  array  in  Tibet.  Three  2D  P- 
wave  velocity  slices  for  the  events  were  used  in  the  SEM  calculations.  All  the  slices  are 
earth- flattened  for  use  in  the  2D  SEM  calculations,  and  the  SEM  results  are  then 
corrected  for  out-of-plane  spreading. 

The  SEM  modeling  of  the  pulse  amplitud  es  and  frequencies  suggests  a  variab  le 
upper  mantle  structure  in  which  the  Qiangta  ng  terrane  has  higher  attenuation  than  the 
Lhasa  terrane  in  the  upper  m  antle.  However,  the  modeling  with  a  constant  upper  m  antle 
attenuation  structure  also  shows  that  th  is  can  provide  a  good  fit  to  the  observed 
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amplitudes  and  frequencies,  but  only  if  a  variable  velocity  gradient  structure  exists  with  a 
lower  velocity  gradient  beneath  the  Qiangtang  terrane  compared  to  the  Lhasa  terrane. 

From  the  observations  and  m  odeling  of  the  seismic  attributes  in  this  study,  two 
alternate  models  for  the  regi  on  are  found,  one  with  a  cons  tant  upper  m  antle  velocity 
gradient  similar  to  those  of  Gr  iffin  et  al.  (2011)  and  Phillips  et  al.  (2007)  with  a  variable 
attenuation  model  and  one  with  a  constant  attenuation  model  and  a  variable  upper  mantle 
velocity  gradient  similar  to  that  of  Myers  et  al.  (2010).  In  either  case,  the  upper  mantle  of 
the  Qiangtang  terrane  can  be  distingui  shed  from  the  Lhasa  terrane.  The  lower  Pn 
velocities,  and  lower  Pn  amplitudes,  suggest  hotter  tern  peratures  in  the  upper  m  antle 
beneath  the  Qiangtang  terrane  in  the  north  co  mpared  to  the  Lhasa  terrane.  Either  lower 
Qp  values  or  a  lower  upper  mantle  velocity  gradient  in  the  Qiangtang  terrane  compared  to 
the  Lhasa  terrane  can  be  used  to  match  the  seismic  attribute  data  presented  here. 


Longitude  (degrees) 

Figure  1.  Map  of  the  stations  of  the  Hi-CLIMB  array  in  Tibet  (triangles)  and  the  epicenters 
of  the  selected  regional  earthquake  sources  (circles).  Black  circles  are  the 
hypocenters  reports  by  the  PDE  catalogs  and  Engdahl  et  al.  (1998),  and  the  white 
circles  are  the  epicenters  nearest  to  the  array  that  have  been  relocated  in  this 
study.  The  dashed  lines  indicate  the  locations  of  north-south  slices  of  the  P-wave 
velocities  from  travel-time  modeling,  with  the  East  Line  at  84  degrees,  the  Central 
Line  at  85  degrees,  and  the  West  Line  at  86  degrees  in  longitude. 

3.  TRAVEL-TIME  MODELING  OF  REGIONAL  EVENTS  IN  TIBET 

Regional  events  recorded  by  the  Tibe  t  Hi-CLIMB  array  with  locations 
approximately  in-line  with  the  array  have  been  chosen  for  analysis  here  and  are  shown  in 
Figure  1 .  However,  other  events  in  specific  te  rranes  have  also  been  id  entitled  for  later 
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processing.  Four  of  the  in-line  events  have  epicenters  less  than  200  km  from  the  Hi- 
CLIMB  array,  and  these  events  were  relocated  with  data  from  the  Hi-CLIMB  array  using 
the  earthquake  location  program  HYPOINVERSE-2000  (Klein,  2002).  The  Engdahl  et  al. 
(1998)  or  the  Preliminary  Determination  of  Epicenters  (PDE)  catalogs  were  used  for  all 
event  locations  excep  t  the  four  nearby  events  (L2a,  L2d,  L3a,  and  L3b),  which  were 
relocated  using  Hi-CLIMB  data. 
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Figure  2.  Ray  diagrams  calculated  using  the  3D  complete  ray  tracing  code  (CRT)  for  the 
Event  LIE  located  by  the  dot. 

The  relocations  improved  the  fit  to  the  travel-times  but  moved  the  epicenters  on  average 
about  5  km.  For  more  distant  events  in  Tibet,  a  focal  depth  normalization  procedure  was 
used  since  events  on  the  Tibetan  plateau  are  dominantly  less  than  15  km  in  focal  depth 
(e.g.,  Molnar  and  Chen,  1983;  Chen  and  Molnar,  1983;  Langin  et  al.,  2003). 

Figure  2  shows  ray  diagram  s  calculated  usi  ng  the  3D  com  plete  ray  tracing  code  CRT 
(Cerveny  et  al.,  1988)  for  the  Event  Lie  shown  by  the  small  circle.  The  left  plot  is  viewed 
from  the  west,  and  the  right  plot  is  v  iewed  from  the  south.  The  blue  top  surface  indicates 
the  structure  of  the  Moho  used  for  the  3D  model.  The  top  plot  of  Figure  3  shows  the 
observed  P-wave  travel  times  for  the  relocate  d  nearby  events  L2a,  L2d,  L3a,  and  L3b, 
compared  with  the  calculated  travel  times  using  an  average  ID  crustal  velocity  m  odel. 
The  bottom  plot  of  Figure  3  shows  the  observed  travel  times  of  the  first-arrival  P  waves 
for  the  nearby  events,  com  pared  with  the  calcu  lated  travel  times  of  a  r  efined,  laterally 
varying  crustal  velocity  model.  Figure  4  shows  the  observed  P-wave  travel-times  (circles) 
for  the  far-regional  Events  (A)  Lie,  (B)  Lla,  (C)  Lli,  (D)  Llj,  (E)  Llk,  and  (F)  E3b.  The 
Pn  travel-times  for  these  events  were  used  to  constrain  the  Moho  structure  and  P-wave 
velocities  in  the  upper  mantle  for  this  part  of  the  Tibetan  Plateau.  The  lines  and  crosses  in 
each  of  the  far-regional  travel-time  plots  are  the  calculated  travel-times  using  the  inferred 
3D  regional  velocity  model. 
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Travel-Times  for  Nearby  Events  L2a,  L2d,  L3a,  and  L3b 
Average  ID  Crustal  Velocity  Model 


Travel-Times  for  Nearby  Events  L2a,  L2d,  L3a,  and  Lib 
Refined  Crustal  Velocity  Model 


Figure  3.  These  plots  show  the  observed  first-arrival  travel  times  and  the  calculated  travel 
times  (lines  and  crosses)  for  the  relocated  nearby  events  using  the  average  ID 
crustal  model  on  top  and  the  refined  crustal  velocity  model  on  the  bottom. 

Figure  5  shows  plots  of  the  3D  Moho  a  nd  upper-mantle  velocity  structure  found  in 
this  study  from  the  modeling  of  near  and  far  re  gional  travel-time  data  for  this  region  of 
the  Tibetan  Plateau,  as  viewed  from  (A)  the  west  and  (B)  the  south  (Griffin  et  al.,  2011). 
The  upper-mantle  velocities  are  higher  and  the  Moho  deeper  beneath  the  southern  portion 
of  the  plateau  south  of  the  Bangong-Nuijang  su  ture  (BNS),  and  upper-mantle  velocities 
are  lower  and  the  Moho  is  shallower  north  of  the  BNS.  North- south  cross-sections  of  the 
3D  velocity  model  of  Griffin  et  al.  (2011)  are  shown  in  Figure  6  for  longitudes  of  84,  85 
and  86  degrees. 

The  structural  results  from  Griffin  et  al.  (201 1)  based  on  the  modeling  of  the  travel¬ 
time  modeling  of  regional  events  are  sim  ilar  to  those  found  by  Tseng  et  al.  (2009)  using 
virtual  seismic  profiles  from  SsPmP  teleseismic  waves  and  by  Nowack  et  al.  (2010)  using 
Gaussian  beam  migration  of  teleseis mic  receiver  functions.  Other  stru  cture  controls  on 
the  structure  have  been  found  by  Hung  et  al.  (2010)  using  finite  frequency  tomography  of 
teleseismic  arrival  times,  and  for  anisotropic  structure  of  the  upper  mantle  by  Chen  et  al. 
(2010). 
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Figure  4.  These  plots  shows  observed  P-wave  travel  times  (circles)  and  calculated  travel 
times  (lines  and  crosses)  for  the  far-regional  events,  (A)  Lie,  (B)  Lla,  (C)  Lli,  (D) 
Llj,  (E)  Llk,  and  (F)  E3b.  The  horizontal  bars  highlight  the  locations  of  the 
transition  for  shoaling  of  the  Moho. 


Figure  5.  The  3D  Moho  and  upper-mantle  velocity  structure  found  in  this  study  for  the 
region  of  the  Hi-CLIMB  array,  with  the  left  plot  showing  a  view  from  the  west  and 
the  right  plot  showing  a  view  from  the  south. 
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Figure  6.  North-south  cross-sections  of  the  3D  velocity  model  of  Griffin  et  al.  (2011)  along 
longitudes  of  84, 85  and  86  degrees. 


4.  SPECTRAL  ELEMENT  MODELING  OF  THE  PN  PHASE 

The  Pn  phase  is  typically  not  a  simple  head  wave  but  is  composed  of  a  c  omplex  set 
of  interfering  waves  when  the  re  is  a  positiv  e  velocity  gr  adient  with  r  espect  to  d  epth 
below  the  Moho  interface  (Cerveny  and  Ra  vindra,  1971;  Hill,  19  71;  Menke  and 
Richards,  1980).  These  interfering  waves  can  also  simply  result  from  the  spherical  nature 
of  the  Earth  (Hill,  1973 ;  Sereno  and  Given,  1990,  Yang  et  al.,  2007),  which  gives  an 
effect  equivalent  to  a  velocity  gradient  in  a  corresponding  flatte  ned  Earth  model.  How 
observed  amplitudes  of  refracted  waves  change  over  distance  supports  the  interpretation 
that  these  arrivals  are  composed  of  diving  and  interference  waves  and  is  in  contrast  to  a 
pure  head  wave  (Braile  and  Sm  ith,  1975).  To  account  for  the  intricate  natu  re  of  Pn 
propagation,  we  have  implemented  the  spectral  element  method  (SEM)  (Komatitsch  and 
Vilotte,  1998;  Komatitsch  etal.,  2005).  SEM  canb  e  applied  in  2D  or  3D  and 
implemented  in  pa  rallel  using  message  passing  interface  (MPI)  on  large-cluster 
computing.  Here  we  have  applied  a  2D  viscoe  lastic  SEM  code.  For  validity  tests  of  the 
SEM  applied  to  the  Pn  phase,  see  Appendix  B.  As  a  practi  cal  measure,  essential  features 
of  Pn  waveforms  are  distilled  into  seismic  attributes,  including  arrival  times,  envelopes  of 
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wave  amplitudes,  and  pulse  frequencies  for  modeling  using  the  approach  of  Matheney 
and  Nowack  (1995)  (see  Appendix  A). 


P-wave  Velocity  Flat  Model  for  Gradient  .004 


Distance  (km) 


P-wave  Velocity  from  Travel-Time  Modeling 


Distance  (km) 


Figure  7.  The  top  plot  shows  the  P-wave  velocities  for  a  simple  layer  over  a  gradient  with  an 
upper-mantle  velocity  gradient  of  .004  1/s.  The  bottom  plot  shows  a  simplified  P- 
wave  velocity  model  from  travel-time  modeling  of  regional  events  from  Hi- 
CLIMB,  showing  a  step  in  the  Moho. 

To  test  the  effects  of  upper-m  antle  velocity  gradient,  Moho  topography,  and  lateral 
velocity  variation  on  Pn  propagation,  we  designed  two  velo  city  models.  The  top  plot  of 
Figure  7  shows  the  P-wave  velocity  m  odel  for  a  laterally  homogeneous  flat  model  with 
an  upper-mantle  velocity  gradient  of  .004  1/s  and  has  a  crustal  layer  above  a  half  space  , 
with  the  Moho  at  70  km  depth.  The  bottom  plot  of  Figure  7  shows  a  simplified  2D  slice 
of  a  model  derived  from  the  3D  travel-tim  e  modeling,  using  regional  data  from  the  Hi- 
CLIMB  array.  Both  the  crust  and  the  upper  mantle  have  velocity  gradients  with  depth  and 
velocities  that  change  laterally.  The  Moho  changes  in  depth  from  73.4  km  on  the  left  side 
of  the  model  to  64.3  km  on  the  right  side  of  the  model. 
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SEM  Synthetic  Travel  Times  for  Different  Gradient  Flat  Models 


SEM  Synthetic  Amplitudes  for  Different  Gradient  Flat  Models 


SEM  Cent  Freq  for  Different  Gradient  Flat  Models 


Distance  (km) 


Figure  8.  SEM  synthetic  P-wave  A)  travel-times,  B)  amplitudes  and  C)  pulse  centroid 
frequencies  for  a  flat  Moho  with  several  velocity  models  with  different  upper- 
mantle  velocity  gradients  of  0.0,  0.002,  0.004,  and  0.006  1/s.  For  this  case  no  Q  is 
included. 

Figure  8  shows  the  first  arrival  amplitudes  for  several  laterally  homogeneous  models 
with  different  upper-m  antle  velocity  gradients  of  0,  0.  002,  0.004,  and  0.006  1/s.  The 
incident  pulse  is  a  Ricker  wavelet  with  a  center  frequency  of  2  Hz  and  filtered  from  .5  to 
2.5  Hz.  The  case  where  the  upp  er-mantle  gradient  is  .004  1/s  is  the  model  shown  in 
Figure  7  (upper  plot).  As  seen  in  Figure  8,  when  there  is  no  upper  m  antle  velocity 
gradient,  the  Pn  pulse  amplitudes  decrease  with  di  stance  from  300  to  800  km  .  For  the 
existence  of  an  upper  mantle  gradient,  the  travel-time  curves  now  curve,  but  the  effect  is 
the  most  pronounced  for  the  highest  gradients.  For  an  upper-m  antle  velocity  gradient 
greater  than  .004,  the  Pn  pulse  amplitudes  increases  for  this  distance  range.  The  pulse 
centroid  frequencies  go  down  at  the  intersection  of  the  Pg  and  Pn  branches,  but  then 
recover  and  even  oscillate  with  distance  for  the  highest  upper  mantle  velocity  gradients. 
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SEM  Synthetic  Travel  Times  for  Different  Gradient  Flat  Models 


SEM  Synthetic  Amplitudes  for  Different  Gradient  Flat  Models 


SEM  Cent  Freq  for  Different  Gradient  Flat  Models 


Figure  9.  SEM  synthetic  P-wave  A)  travel-times,  B)  amplitudes  and  C)  pulse  centroid 
frequencies  for  a  flat  Moho  with  several  velocity  models  with  different  upper- 
mantle  velocity  gradients  of  0.0,  0.002,  0.004,  and  0.006  1/s.  For  this  case  a  Q  of  150 
is  assumed  in  the  crust  and  400  in  the  mantle. 

In  Figure  9,  a  similar  case  as  in  Figure  8  is  shown  with  a  Q  of  150  in  the  crust  and 
400  in  the  m  antle.  The  travel-tim  es  are  th  e  same  as  befor  e  and  still  bend  down  for 
increasing  upper  mantle  velocity  gradient.  Ho  wever,  the  attenuation  has  the  effect  of 
bringing  down  all  the  am  plitude  curves.  None  theless,  the  amplitude  pattern  is  still 
preserved.  The  pulse  centroid  frequencies  are  also  brought  down. 

Figure  10  shows  the  travel-tim  es,  amplitudes  and  pulse  centroid  f  requencies 
computed  for  the  velocity  m  odel  shown  in  the  lower  plot  of  Figure  7,  with  a  step  in  the 
Moho.  The  crosses  are  for  sour  ce  B  near  1000  km  and  recorded  to  the  south  in  Figure  7 
(right  plot),  and  the  circles  are  for  source  A  near  0  km  and  recorded  to  the  north  in  Figure 
7  (right  plot).  The  trav  el-time  curves  are  qui  te  distinct  be  tween  the  reversed  events 
resulting  from  the  different  upper  mantle  top  velocities  to  the  left  and  right  of  the  model, 
as  well  as  the  change  of  Moho  structure  acros  s  the  m  odel.  The  am  plitudes  curves  in 
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Figure  10  show  a  relative  increase  with  distance  for  the  Pn  branch  for  source  B,  resulting 
from  the  step  in  the  Moho  for  this  particul  ar  upper-mantle  velocity  gradient.  However, 
for  source  A,  the  Pn  amplitudes  relatively  decrease  with  distance,  resulting  from  the  step 
in  the  Moho.  However,  the  effects  of  the  Moho  step  on  the  pulse  frequencies  are  less 
dramatic.  In  any  case,  using  events  at  different  locations  to  the  north  and  south  of  the  Hi- 
CLIMB  array,  one  can  identify  the  large-scale  structural  features. 


Travel  Times  for  Model- 1 1  Normal  and  Reverse  Shots 
Source  =  Isotropic  Moment  Tensor,  f0=  1.75  Hz 


Amplitudes  vs  Distance  for  Model- 1 1  Normal  and  Reverse  Shots 
Source  =  Isotropic  Moment  Tensor,  to=  1.75  Hz 


Distance  (km) 


Cent  Freq  for  Model- 11  Normal  and  Reverse  Shots 
Source  =  Isotropic  Moment  Tensor,  f0  =  1 .75  Hz 


Figure  10.  SEM  synthetic  A)  travel-times,  B)  amplitudes  and  C)  pulse  centroid  frequencies 
for  a  step  in  the  Moho,  as  shown  for  the  right  plot  in  Figure  7.  The  crosses  for 
source  B  near  1000  km  in  Figure  7  (right  plot)  and  recorded  to  the  south,  and  the 
circles  are  for  a  source  A  near  0  km  and  recorded  to  the  north  in  Figure  7  (right 
plot). 

Figure  1 1  shows  the  case  of  a  step  in  the  Moho  given  in  the  right  plot  in  Figure  7, 
only  now  with  attenuation  included.  The  crust  now  has  a  Q  of  150  and  the  upper  m  antle 
has  a  Q  of  400.  Once  again  the  travel-times  are  the  same  as  the  no  attenuation  case  shown 
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in  Figure  10.  For  the  attenuation  case  in  Figure  1 1,  the  amplitudes  are  brought  down  but 
with  a  similar  pattern  as  in  the  no  attenu  ation  case.  Also,  the  effect  of  the  step  in  the 
Moho  on  pulse  frequencies  is  m  ore  subdued  in  the  case  with  attenuation.  However,  the 
Pn  pulse  frequencies  still  show  a  slight  increase  with  frequencies  resulting  from  the  upper 
mantle  velocity  gradient. 


Travel  Times  for  Model- 1 1  Normal  and  Reverse  Shots 
Source  =  Isotropic  Moment  Tensor,  f0=  1.75  Hz 


Amplitudes  vs  Distance  for  Model- 1 1  Normal  and  Reverse  Shots 
Source  =  Isotropic  MT,  Q  crust  =  1 50  and  Q  mantle  =  400 


Distance  (km) 


Cent  Freq  for  Model- 1 1  Normal  and  Reverse  Shots 
Source  =  Isotropic  MT,  Q_crust  =  1 50  and  Q_mantlc  =  400 


Figure  11.  SEM  synthetic  A)  travel-times,  B)  amplitudes  and  C)  pulse  centroid  frequencies 
for  a  step  in  the  Moho,  as  shown  for  the  right  plot  in  Figure  7.  The  crosses  for 
source  B  near  1000  km  in  Figure  7  (right  plot)  and  recorded  to  the  south,  and  the 
circles  are  for  a  source  A  near  0  km  and  recorded  to  the  north  in  Figure  7  (right 
plot). 

5.  MOMENT  TENSORS  OF  REGIONAL  EVENTS  IN  TIBET 

We  have  compiled  the  centroid  m  oment  tensor  solutions  of  regional  earthquakes 
that  are  relevant  to  those  selected  for  constructing  Pn  profiles  from  CMT  focal 
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Event  ID:  E3b 

Date:  2004  /  1 2  /  5  Time:  23:35: 38.8 
Epicenter:  34.3700  N  86.2400  E 
Depth:19  km 
Mw:3.817 

Nodal  plane:  236.23  /  71 .1 5  / 14.41 
Nodal  plane:  141.49/76.38/ 160.58  a 
CLVD:  30.64  %  Isotropic:  0.00  %  isk 
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Average  misfit  for  21  stations  =  0.523 


Figure  12.  Inversion  of  three-component  waveforms  for  the  focal-mechanism  of  event  E3b 
using  data  from  the  Hi-CLIMB  array.  The  top  right  plot  shows  the  data  misfit 
with  a  narrow  minimum  at  6  km  and  and  a  broad  flat  minimum  at  19  km,  and  the 
upper  left  plot  shows  the  19  km  solution.  The  lower  plots  show  the  matches 
between  the  observed  and  calculated  waveform  data. 


mechanisms  of  Dziewonski  et  al.  (1981)  and  from  Baur  (2007).  However,  event  E3b  and 
some  of  the  smaller  events  did  not  have  previously  derived  focal  mechanisms.  For  the  Hi- 
CLIMB  array-station  geometry  for  some  of  the  regional  events,  making  a  determination 
of  the  focal  mechanisms  is  relatively  challenging  even  though  the  data  waveform  quality 
is  very  high  from  the  Hi-CLIMB  array.  As  such,  Tai-Lin  Tseng  of  our  group  checked 
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with  the  regional  LS A  s  tation  (the  only  perm  anent  regional  seismic  station)  and  other 
temporary  arrays  in  Tib  et  for  their  availabil  ity  and  tried  to  include  th  em  in  the  m  oment 
tensor  inversions.  She  found  that  for  a  number  of  earthquakes  investigated  they  are  either 
shallow  normal  faults  on  the  rifts  or  strike  -slip  faults  occurring  in  the  upper  crust. 
However,  both  types  of  mechanisms  show  generally  E-W  extensional  stresses  which  are 
typical  for  the  Tibet  area. 

As  an  example  of  a  focal  mechanism  inversion,  Figure  12  shows  the  determ  ination 
of  the  focal  mechanism  for  event  E3b  using  seismic  waveform  data  recorded  by  the  Hi- 
CLIMB  array.  The  top  right  plot  shows  the  data  misfit  with  a  narrow  m  inimum  at  5  km 
and  a  broad  flat  minimum  at  19  km.  Figure  12  also  shows  the  waveform  matches  between 
the  observed  and  calculated  long  period  data  for  the  19  km  solution.  F  or  this  event,  the 
seismograms  for  the  permanent  seismic  station  LSA  were  poor. 


Range  of  Focal  Mechanism  Solutions  for  Event  E3b 


Figure  13.  Focal  mechanism  solutions  for  Event  E3b.  a)  A  comparison  global  CMT  focal 
mechanism  for  an  earthquake  located  about  70  km  distant  from  Event  E3b.  b) 
The  focal  mechanism  solution  for  Event  E3b  with  a  data  misfit  minimum  at  a 
depth  of  6  km.  c)  The  focal  mechanism  solution  for  Event  E3b  with  a  data  misfit 
minimum  at  a  depth  of  19  km.  d)  The  focal  mechanism  for  Event  E3b  for  a  shift 
in  epicenter  of  10  km  to  the  east  resulting  in  a  least  misfit  minimum  at  the 
shallower  depth  of  7  km.  All  solutions  indicate  East-West  extension. 

Figure  13  shows  a  summary  of  the  fo  cal  mechanism  solutions  for  Event  E3b.  As  a 
comparison,  Fig.  13a  shows  a  global  CMT  focal  mechanism  of  an  earthquake  located 
about  70  km  distant  from  Event  E3b.  For  E  3b,  there  are  two  local  minimum,  one  narrow 
at  around  6  km  in  depth  (Fig.  13b),  and  one  br  oad  and  flat  at  around  19  km  depth  (Fig. 
13c).  However,  the  average  m  isfits  for  this  event  are  still  r  datively  large  at  any  depth, 
and  both  solutions  have  high  CLVD  components.  This  initially  suggested  to  us  that  there 
could  be  either  variations  in  the  local  veloc  ity  model  or  event  location  issues  that  c  ould 
need  to  be  accounted  for. 
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An  additional  experiment  was  don  e  to  eh  eck  the  s  ensitivity  to  the  epicentral 
location  and  the  moment  tensor  solutions  were  mostly  stable,  except  for  a  epicentral  shift 
of  about  10  km  to  the  east  which  resulted  in  the  shallow  solution  at  7  km  becoming  the 
solution  with  the  sm  allest  data  misfit  (Fig.  13d),  and  is  similar  to  the  ear  lier  misfit 
minimum  solution  at  6  km.  Nonetheless,  Ev  ent  E3b  has  a  large  az  imuthal  gap  making  it 
a  challenging  event  to  get  a  reliable  focal  mechanism  solution  for.  However  all  solutions 
for  this  event  indicate  o  verall  east-west  extension.  We  are  currently  using  both  inverted 
mechanisms  for  event  E3b  for  the  modeling  of  Pn  seismic  attributes. 

6.  SPECIFICATION  OF  THE  MOMENT  TENSOR  FOR  2D  SEM 
MODELING 

For  the  specification  of  the  moment  tensor  in  the  2D  SE  M  modeling,  we  first  tested 
several  2D  azimuth  slices  of  a  3D  moment  tensor  along  specific  azimuths  for  a  particular 
double-couple  moment  tensor  (left  plot  of  the  figure  14).  The  Pg  and  Pn  take-off  angles 


Radiation  Patterns  for  Different  Azimuthal  Angles 


Take-off  Angle  (Degree) 


Figure  14.  The  left  plot  shows  the  3D  focal  mechanism  for  a  particular  double-couple 
moment  tensor  source  with  four  angles  given  by  dashed  lines.  The  Pg  and  Pn  take¬ 
off  angles  are  shown  by  small  circles.  The  right  plot  shows  the  amplitudes  on  the 
focal  sphere  of  the  2D  moment  tensor  for  the  azimuthal  angles  given  in  the  left  plot 
as  a  function  of  take-off  angle  derived  from  the  3D  moment  tensor. 

are  shown  by  small  circles.  The  radiation  pattern  terms  resulting  from  implementing  a 
rotated  2D  mom  ent  tensor  derived  from  the  3D  moment  tensor  are  then  shown  in  the 
right  plot  of  figure  14. 

We  next  com  pare  the  P-wave  amplitude  s  with  the  amplitudes  from  the  2D 
Specfem2D  code  using  the  rotated  2D  m  oment  tensor  formulas.  Based  on  these  tests,  we 
inferred  that  an  additional  180  0  azimuthal  rotation  is  needed  to  m  atch  the  amplitudes  of 
the  rotated  moment  tensor.  With  this  modification,  the  Pn  and  Pg  amplitude  computed 
from  the  SEM  code  are  shown  in  figure  1 5  for  the  focal  mechanisms  given  the  plot  to  the 
left  using  four  azimuthal  angles. 
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7.  MODELING  PN  SEISMIC  ATTRIBUTES  IN  TIBET 


Numerical  modeling  of  Pn  wave  attributes,  including  the  travel-times,  the  peak 
envelope  amplitudes,  and  the  instantaneous  pulse  frequencies,  was  performed  for  six 


SEM  Envelope  Amplitudes  for  Different  Angles 


□ 


Figure  15.  The  left  plot  shows  a  focal  mechanism  with  four  azimuthal  angles  given  by 
dashed  lines.  The  Pn  and  Pg  take-off  angles  are  shown  by  the  small  circles.  The 
right  plot  shows  the  SEM  envelope  amplitudes  for  the  azimuthal  angles  given  on 
the  left.  Different  symbols  represent  different  azimuthal  angles. 


Focal  Mechanism  for  L I E  with  Hi-Climb  Stations  Focal  Mechanism  for  E3B  with  Hi-Climb  Stations  Focal  Mechanisms  for  South  Events  with  Hi-Climb  Stations 


Longitude  (Degree)  Longitude  (Degree)  Longitude  (Degree) 

Figure  16.  The  left  map  shows  the  original  and  slightly  revised  focal  mechanism  for  event 
LIE.  The  middle  map  shows  the  focal  mechanism  TL-1  and  an  alternative  focal 
mechanism  TL-2  for  event  E3B.  The  right  map  shows  the  focal  mechanisms  for  the 
cluster  of  south  events  L1A,  L1I,  L1J,  and  L1K.  The  triangles  show  the  Hi-Climb 
stations,  where  the  solid  triangles  along  with  the  solid  lines  show  the  orientation 
angles  used  for  the  specification  of  the  moment  tensor.  The  dashed  lines  show  the 
separation  of  ranges  used  to  specify  the  moment  tensor. 
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regional  earthquakes  recorded  by  the  Hi-Climb  array  in  Tibet.  Earthquakes  are  chosen  to 
be  approximately  in-line  with  the  Hi-Climb  array,  but  also  include  both  Pg  and  Pn 
arrivals  along  with  the  cross-over  distance.  The  location  o  f  the  events  used  here  are 
shown  in  figure  16. 


Figure  17.  This  map  shows  the  locations  where  2D  P-wave  velocity  slices  were  obtained 
from  the  3D  velocity  model  of  Griffin  et  al.  (2011).  Triangles  represent  the  Hi-Climb 
array.  The  circles  show  the  events  used  by  Griffin  et  al.  (2011)  for  the  3D  velocity 
model  where  the  open  circles  are  relocations.  The  dashed  lines  show  the  West, 
Central,  and  East  lines  for  the  2D  slices. 

For  the  SEM  synth  etic  calculations  using  the  selected  Hi-Climb  events,  P-wave 
velocities,  upper  mantle  velocity  gradients,  and  the  Moho  structures  from  Griffin  et  al. 
(201 1)  was  used.  The  Moho  structure  found  by  Griffin  et  al.  (2011)  using  regional  travel- 
times  is  similar  to  that  found  by  Tseng  et  al.  (2009)  and  Nowack  et  al.  (2010)  using 
teleseismic  waves.  The  upper  mantle  top  velocities  and  velocity  gradients  from  Griffin  et 
al.  (2011)  are  sim  ilar  to  those  found  by  Phillip  s  et  al.  (200  7).  These  are  also  sim  ilar  to 
those  found  by  Myers  et  al.  (2010)  in  the  so  uthem  Lhasa  terrane,  but  the  velocity 
gradients  from  Myers  et  al.  (2010)  are  lower  in  the  northern  Qiangtang  terrane. 

Griffin  et  al.  (2011)  used  the  3D  ray  tr  acing  code  CRT  of  Cerveny  et  al.  (1988) 
for  modeling  the  travel-times  and  imaged  a  3D  P-wave  velocity  m  odel  along  with  the 
Moho  structure.  From  this  3D  model,  2D  slices  were  obtained  for  use  in  the  2D  spectral- 
element  method  (SEM)  m  odeling.  The  location  of  the  2D  slices  and  the  3D  P-wave 
velocity  model  of  Griffin  et  al.  (201 1)  ar  e  given  in  figure  17.  On  e  2D  slice  along  the 
longitude  of  84°  (West  Line)  is  co  nsidered  for  the  south  event  LIE  for  the  seis  mic 
attribute  modeling.  Another  2D  slice  along  the  longitude  of  85  0  (Central  Lin  e)  is 

considered  for  the  south  cluster  of  events  LI  A,  LI  I,  LI  J,  an  d  L1K.  The  2D  slice  along 
the  longitude  of  86°  (East  Line)  is  c  onsidered  for  the  north  event  E3B.  Figure  18  shows 
the  P-wave  velocity  s  lices  for  the  event  LIE,  the  event  E3  B,  and  the  cluster  of  events 
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LI  A,  LI  I,  LI  J,  and  L1K  .  These  velocity  pr  ofiles  are  then  earth  flattened  for  input  into 
the  SEM  code  (see  Bakir,  2011). 


P-Wave  Velocity  Model  for  Event  LIE 


P-Wave  Velocity  Model  for  Event  E3B 

Distance  (km)  N 

200  400  Ann  800  1000 


P-Wave  Velocity  Model  for  Event  LI  A,  LI  I,  L1J,  and  L1K 
s  Distance  (km)  N 


Figure  18.  The  plots  show  the  P-velocities  for  3  different  slices  of  the  3D  model  of  Griffin  et 
al.  (2011)  which  have  constant  upper  mantle  velocity  gradients  but  different  upper 
mantle  top  velocities.  The  small  circles  show  the  locations  of  the  events  used  with 
each  velocity  slice. 

The  moment  tensor  solutions  of  the  far  regional  6  events  were  obtained  from  Baur 
(2007)  for  all  events  except  E3B  which  was  obtained  from  Tai-Lin  Tseng  from  our  group 
and  described  earlier.  For  the  numerical  calculations,  SPECFEM2D,  a  parallel  SEM  code 
(Komatitsch  andVilotte,  1998  and  Komatitsch  etal.,  200  5),  was  used  sim  ilar  to  the 
previous  numerical  tests.  For  each  event,  m  ultiple  2D  mom  ent  tensor  sources  are 
calculated  for  different  azimuth  angles  from  the  events  to  the  Hi-C  limb  stations.  Two, 
three,  or  four  stations  are  selected  for  the  orientation  of  the  moment  tensor  depending  on 
how  much  the  focal  m  echanism  solutions  are  affected  by  the  azim  uth  to  the  statio  ns. 
Three  stations  are  chos  en  for  LIE  for  the  specification  of  the  moment  tensor.  Four 
stations  are  chosen  for  E3B  since  it  is  the  most  out-of-line  event  and  stations  get  closer  to 
the  nodal  planes  for  the  different  azimuth  angl  es.  Two  stations  are  found  to  be  sufficient 
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for  the  cluster  of  events  L1A,  L1I,  LI  J,  and  L1K  since  they  are  the  most  in-line  with  the 
Hi-CLIMB  array. 


Distance  (km)  Distance  (km) 


Event  L 1 K  compared  with  No  Q  Event  E3B  TL- 1  compared  with  No  Q 


Distance  (km) 


Distance  (km) 


Figure  19.  These  plots  show  the  calculated  travel-times  of  the  selected  events  with  no 
attenuation  in  the  SEM  models  compared  with  the  observed  travel-times.  Circles 
are  the  SEM  calculated  travel-times  and  crosses  are  observed  travel-times. 

In  all  peak  envelope  amplitude  and  pulse  frequency  figures,  vertical  dashed  lines  are  used 
to  represent  different  distance  ranges  defined  by  the  different  orientation  angles  used  for 
the  2D  moment  tensors.  The  approximate  orientation  angles  and  the  strike,  rake,  and  dip 
of  the  events  are  used  to  calculate  the  m  oment  tensor  entries  of  the  focal  mechanisms. 
Once,  a  rotated  3D  moment  tensor  matrix  is  obtained  as  described  by  Bakir  (2011),  only 
the  x  and  z  com  ponents  of  the  rotated  m  oment  tensor  are  input  into  the  2D  SEM 
calculations.  Figure  16  shows  the  moment  tensors  of  the  selected  events,  the  Hi-CLIMB 
array,  and  the  orientation  angl  es  from  the  ev  ents  to  th  e  selected  stations  on  a  regional 
map  for  the  specification  of  the  2D  moment  tensors. 
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Event  LIE  Envelope  Amplitudes  with  No  Q 
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Figure  20.  These  plots  show  the  peak  envelope  amplitudes  of  the  selected  events  with  no 
attenuation  in  the  SEM  models  compared  with  the  observed  envelope  amplitudes. 
The  circles  are  the  SEM  calculated  values  and  the  crosses  are  observed  envelope 
amplitudes.  The  distance  ranges  for  the  specified  orientation  angles  for  the  moment 
tensor  are  also  shown.  For  event  LIE,  the  modified  moment  tensor  is  used  and  for 
event  E3B  the  TL-1  moment  tensor  is  shown. 


Travel-times  from  the  SEM  elastic  calcula  tions  are  f  irst  matched  with  the  obse  rved 
travel-times  from  regional  earthquakes  reco  rded  by  the  Hi-CLIMB  array  in  figure  19. 
From  figure  19,  it  can  seen  that  the  2D  slices  derived  from  the  3D  model  of  Griffin  et  al. 
(201 1)  do  a  reasonable  job  of  m  atching  the  observed  travel-times.  The  peak  envelope 
amplitudes  and  the  pulse  frequencies  of  the  SEM  calculations  are  then  compared  with  the 
observed  data.  Because  the  SEM  results  will  be  compared  with  the  observed  data,  the  2D 
SEM  amplitudes  are  corrected  for  out-of-p lan  e  spreading  (see  Bakir,  2011).  Both  the 
centroid  pulse  frequencies  and  the  instantane  ous  frequencies  are  used  to  describe  the 
frequency  characteristics  of  the  data  usi  ng  the  approach  described  by  Matheney  and 
Nowack  (1995). 
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Figure  21.  These  plot  show  the  calculated  instantaneous  pulse  frequencies  of  the  selected 
events  with  no  attenuation  for  the  SEM  models  compared  with  the  observed 
instantaneous  and  pulse  frequencies.  The  circles  are  the  SEM  calculations  and  the 
crosses  and  squares  are  the  observed  instantaneous  and  centroid  frequencies.  For 
event  LIE,  the  modified  moment  tensor  is  used  and  for  event  E3B  the  TL-1  moment 
tensor  is  shown. 


Figure  20  shows  the  calculated  peak  e  nvelope  amplitudes  of  the  elas  tic  SEM 
calculations  assuming  no  atte  nuation  compared  with  the  observed  peak  envelope 
amplitudes  for  the  six  events.  For  event  LIE,  the  revised  moment  tensor  is  used  to  avoid 
small  amplitudes  of  the  Pg  branch  being  too  close  the  nodal  plane  for  the  original 
moment  tensor  of  Baur  (2007).  However,  this  was  within  the  error  of  the  determination  of 
the  moment  tensor  (J.  Nabelek,  personal  communication).  For  event  E3B,  there  were  two 
allowed  focal  m  echanisms  as  derived  by  Tai-  Lin  Tseng.  Here,  we  give  the  results  for 
focal  mechanism  TL-1  shown  in  Figure  16.  Ho  wever,  the  results  for  TL-2  were  similar 
and  are  given  by  Bakir  (201 1). 
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Figure  22.  These  plots  show  variable  upper  mantle  Q  models  for  each  slice  determined  from 
modeling  seismic  attributes  from  the  selected  regional  events  assuming  the  velocity 
slices  in  Figure  18. 
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Figure  21  shows  the  pulse  fr  equencies  using  the  elastic  SEM  ca  lculations  with  no 
attenuation  compared  with  the  observed  inst  antaneous  pulse  frequencies.  The  calculated 
pulse  frequencies  for  the  Pn  branch  increase  with  dist  ance  resulting  from  the  upper 
mantle  velocity  gradient  obtained  from  the  3D  models  Griffin  et  al.  (2011)  and  similar  to 
those  obtained  by  Phillips  et  al.  (2007).  The  pulse  frequencies  can  be  reduced  by  either 
increasing  the  attenuation  or  lowering  the  velocity  gradient  in  the  upper  m  antle. 
However,  both  th  e  calculated  envelope  amplitudes  and  pulse  f  requencies  are  too  h  igh 
when  using  no  attenuation  for  these  models  comp  ared  to  the  observed  data.  From  this  we 
conclude  that  attenuation  or  a  lower  gradient  in  the  upper  mantle  is  needed  to  m  atch  the 
amplitudes  and  pulse  frequencies  with  the  observed  attributes. 
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Figure  23.  These  plot  show  the  SEM  calculated  peak  envelope  amplitudes  of  the  selected 
events  with  the  attenuation  models  shown  in  Figure  22  and  the  velocity  slices  given 
in  Figure  18  compared  with  the  observed  envelope  amplitudes.  The  circles  are  the 
SEM  calculated  values  and  the  crosses  are  observed  envelope  amplitudes.  The 
distance  ranges  for  the  specified  orientation  angles  for  the  moment  tensor  are  also 
shown.  For  event  LIE,  the  modified  moment  tensor  is  used  and  for  event  E3B  the 
TL-1  moment  tensor  is  shown. 

From  the  modeling  of  the  am  plitudes  and  pulse  frequencies  given  in  figures  20, 
and  21,  this  suggests  that  uppe  r  mantle  attenuation  is  approp  riate  to  apply  for  all  the 
events.  For  the  modeling  here,  we  first  use  the  structural  models  obtained  from  Griffin  et 
al.  (2011)  and  then  include  seismic  attenuation.  A  constant  Q  attenuation  with  frequency 
is  provided  in  the  spectral  elem  ent  code,  and  given  that  the  data  and  the  synthetics  are 
both  filtered  in  the  band  between  .5  to  2.5  Hz  is  probably  ap  propriate.  However,  further 
tests  of  this  will  be  required  in  the  future. 
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Figure  24.  These  plot  show  the  SEM  calculated  pulse  frequencies  of  the  selected  events  with 
the  attenuation  models  shown  in  Figure  22  and  the  velocity  slices  given  in  Figure  18 
compared  with  the  observed  pulse  frqeuencies.  The  circles  are  the  SEM  calculations 
and  the  crosses  and  squares  are  the  observed  instantaneous  and  centroid 
frequencies.  For  event  LIE,  the  modified  moment  tensor  is  used  and  for  event  E3B 
the  TL-1  moment  tensor  is  shown. 

For  the  southwest  event  LIE,  a  Qp  of  130  in  the  crust,  a  Q  p  of  260  in  the  Lhasa 
terrane  to  the  south  and  a  Qp  of  130  in  the  Qiangt  ang  terrane  in  the  north  for  the  upper 
mantle  was  used  for  the  SEM  m  odeling.  For  the  south  cluster  even  ts,  Qp  values  were 
inferred  to  be  the  sam  e  as  event  L  IE,  except  a  Qp  of  100  is  used  in  the  crust.  For  the 
north  event  E3B,  a  Qp  of  100  in  the  crust  and  a  Qp  of  260  for  the  Lhasa  terrane  in  the 
south  and  a  Qp  of  120  for  the  Qiangtang  terr  ane  in  the  north  was  used  in  the  upper 
mantle.  When  a  variable  Qp  exists,  a  s  mooth  transition  is  created  around  the  ramp 
structure  of  the  Moho  for  the  SEM  models.  The  variable  Qp  models  of  these  events  are 
shown  in  figure  22. 
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Figure  25.  The  plots  show  the  P-velocities  for  3  different  slices  of  the  3D  model  of  Griffin  et 
al.  (2011)  which  variable  upper  mantle  velocity  gradients.  The  small  circles  show 
location  of  the  events  used  with  each  velocity  slice. 

In  figure  23,  the  SEM  envelope  amplitude  with  variable  attenuation  and  a  constant 
velocity  gradient  in  the  upper  m  antle  are  compared  with  the  observed  peak  envelope 
amplitudes  and  in  figure  24  the  com  puted  pulse  frequencies  are  com  pared  with  the 
observed  instantaneous  and  centroid  frequencies.  From  figures  23  and  24,  it  is  observed 
that  the  SEM  m  odeling  with  variable  attenuation  and  a  cons  tant  upper  mantle  velocity 
gradient  gives  a  bette  r  fit  to  the  seismic  attributes  than  the  earlier  case  with  no 
attenuation. 
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Figure  26.  The  plots  shows  constant  upper  mantle  attenuation  used  with  the  velocity  models 
in  Figure  25  with  a  laterally  variable  upper  mantle  velocity  gradient  between  the 
Lhasa  and  Qiangtang  terranes. 

We  next  investigate  if  ve  locity  models  with  variab  le  upper  m  antle  velocity 
gradients  can  give  sim  ilar  results  for  the  m  odeling  of  seismic  attributes  for  the  d  istance 
ranges  used  here.  A  variable  upper  m  antle  gradient  structure  is  constructed  for  all  the  2D 
slices  and  input  into  the  SEM  calculations.  A  constant  Qp  of  260  is  chosen  for  the  entire 
upper  mantle  along  with  a  vari  able  upper  mantle  velocity  gradient.  Figure  25  shows  the 
variable  gradient  P-wave  velocity  slices  for  the  southeast  event  LIE,  the  north  event  E3B, 
and  for  the  south  cluster  events  L  1A,  LI  I,  L1J,  and  L  IK.  W  hile  the  earlier  P-wave 
velocity  slices  had  a  velo  city  gradient  of  0.003  sec  1  over  the  entire  upper  m  antle  before 
the  earth-flattening,  the  variable  velocity  gradient  slices  have  a  velocity  gradient  of  0.003 
sec'1  in  the  Lhasa  terrane  to  the  south  and  a  velocity  gradient  of  0.001  sec'1  in  the 
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Figure  27.  These  plots  show  the  peak  envelope  amplitudes  of  the  selected  events  with  a 
constant  attenuation  in  the  upper  mantle  for  the  SEM  models  compared  with  the 
observed  envelope  amplitudes.  The  circles  are  the  SEM  calculations  and  the  crosses 
are  observe  envelope  amplitudes. 

Qiangtang  terrane  to  th  e  north.  This  is  sim  ilar  to  the  upp  er  mantle  velocity  grad  ient 
model  of  Myers  et  al.  (2010)  for  this  region.  The  corresponding  Q  models  with  constant 
upper  mantle  attenuation  are  shown  in  figure  26. 

For  all  eve  nts,  the  vis  coelastic  SEM  results  with  constant  attenuation  and  a 
variable  upper  mantle  gradient  structure  are  then  com  pared  with  the  observed  data.  In 
figure  27,  a  comparison  of  the  SEM  amplitudes  with  constant  attenuation  and  a  variable 
upper  mantle  gradient  with  the  observed  amplitudes  is  given,  an  d  in  figure  28  a 
comparison  of  the  SEM  pulse  frequencies  wi  th  the  observed  instan  taneous  and  centroid 
pulse  frequencies  is  shown.  The  m  odeling  results  are  similar  to  the  results  obtained  for 
the  variable  attenuation  and  constant  gradient  results.  This  suggests  a  trade-off  in  model 
parameters  between  upper  m  antle  attenuation  and  upper  mantle  velocity  gradients  using 
the  amplitude  and  pulse  frequency  data. 
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Figure  28.  These  plots  show  the  instantaneous  pulse  frequencies  of  the  selected  events  with  a 
constant  attenuation  and  a  variable  upper  mantle  velocity  gradient  for  the  SEM 
models  compared  with  the  observed  instantaneous  and  centroid  frequencies.  The 
circles  are  the  SEM  calculations  and  the  crosses  and  squares  are  observed 
instantaneous  and  centroid  frequencies. 


8.  SEISMIC  PROFILING  WITH  TELESEISMIC  WAVES 

In  an  effort  to  better  constrain  the  structural  model  for  the  re  gion,  we  have  also 
completed  several  studies  of  large-scale,  north-south  variations  in  crustal  thickness  across 
the  Tibetan  plateau  using  teleseism  ic  waves.  Tseng  et  al.  (2009)  devised  a  new  approach 
using  the  SsPmp  phase  for  constructing  deep-pen  etrating  seismic  profiles  f  rom 
teleseismic  data  which  revea  1  significant,  regional  variations  in  crustal  thickness  under 
near-constant  elevation  of  Tibet. 
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Figure  29.  North-south  trending  profile  across  southern  and  central  Tibet,  showing 
topography  on  top,  a  P-wave  reflection  profile  constructed  using  the  phase  SsPmp, 
and  the  estimated  crustal  thickness  (Panels  a,  b,  c,  respectively). 

The  SsPmp  phase  originates  as  the  direct  S'- wave  (Ss  phase)  reflects  under  the  free 
surface  and  partially  converts  to  P-  wave,  which  in  turn  reflects  off  the  top  of  the  Moho 
before  finally  arriving  at  a  seismic  station.  The  key  concept  is  that  the  P-  to  S-  conversion 
point  below  the  free  surface  is  a  virtual  source  for  the  subsequent,  reflected  P- wave.  In 
this  view,  the  SsPmp  phase  is  equivalent  to  the  prominent  “PwP”  phase  in  a  conventional 
seismic  profile.  Considering  all  known  sources  of  error  in  this  approach,  the  overall  final 
uncertainty  of  crustal  thickness  is  no  more  than  about  ±3  km  (over  a  total  thickness  of 
about  70  km,  or  ±4%). 

The  very  thick  crust  under  Tibet  is  by  no  means  uniform,  instead  there  is  an  overall 
trend  of  northward  th  inning.  The  Moho  reach  es  its  deep  est  depth,  up  to  73-77  km, 
beneath  the  southern  Lhasa  terrane  but  shoals  to  around  60-64  km  beneath  the  Qiangtang 
terrrane.  This  regional  vari  ation  in  crustal  th  ickness  has  a  profound  effect  on  the 
propagation  of  the  Pn  phase,  and  has  been  incorporated  into  the  three-dim  ensional  model 
of  the  lithosphere  for  detailed  modeling  of  the  frequency-dependent  propagation  of  the 
Pn  phase.  Figure  29  shows  a  north-south  tren  ding  profile  across  southern  and  central 
Tibet,  showing  topography  on  top,  a  pseudo  P-  wave  reflection  profile  constructed  using 
the  SsPmp  phase,  and  the  estimated  crustal  thickness  (Panels  a,  b,  c,  respectively). 

We  also  co  nducted  Gaussian  beam  migration  of  teleseism  ic  receiver  functions 
(Nowack  et  al.,  2010).  This  approach  alio  ws  for  better  resolution  than  traditional 
piercing  point  imaging  using  receiver  functions.  Figure  30  shows  a  comparison  between 
crustal  thickness  estimated  from  wide-angle  P-wave  reflections  from  modeling  the  SsPmp 
phases  (dots)  and  an  im  age  of  the  Tibetan  lithosphere  obtained  from  Gaussian  beam 
migration  of  direct  P-  to  S-wave  conversi  ons  (from  Tseng  et  al.,  2009;  and  the  Gaussian 
beam  migration  from  Nowack  et  al.,  2010). 
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Figure  30.  A  comparison  between  crustal  thickness  estimated  from  wide-angle  P-wave 
reflections  (dots)  from  the  teleseismic  SsPmp  phase  and  an  image  of  the  Tibetan 
lithosphere  obtained  from  Gaussian  beam  migration  of  direct  P-  to  S-wave 
conversions  (from  Tseng  et  al.,  2009;  and  the  Gaussian  beam  migration  from 
Nowack  et  al.,  2010). 

Thinner  crust  under  northern  Tibet  has  been  taken  as  evidence  for  full  isostasy 
being  achieved  below  the  lithosphere,  deep  in  the  upper  m  antle  where  tem  perature  is 
anomalously  high.  Supporting  evidence  for  such  a  m  echanism  includes  less  efficient 
propagation  of  Sn  waves  just  below  the  Moho,  and  larg  e  birefringence  of  S-waves  (up  to 
2  s  or  more  between  fast-  and  slow-directions  of  polarization)  that  has  been  inferred  as  a 
result  of  east-west  flow  in  the  upp  er  mantle.  While  our  findings  are  generally  consistent 
with  this  interpretation  in  a  broad  sense,  there  are  preliminary  indications  that  variations 
in  properties  of  the  upperm  ost  mantle  are  too  abrupt  to  be  accounted  for  by  therm  al 
effects  alone. 

Finally,  Hung  et  al.  (2010)  and  Hung  et  al.  (201 1)  have  developed  in  detail,  both  the 
theory  and  the  practice  of  an  innovation  th  at  combines  finite-frequency  (“banana- 
doughnut”),  travel-time  tomography  and  a  d  ata-adaptive,  multi-scale  approach.  W  hen 
applied  to  the  teles  eismic  P-  and  5-wave  data,  collected  un  der  Project  Hi-CLIMB,  this 
new  method  yields  high-resolution  results  in  the  uppermost  mantle  where  both  lateral  and 
radial  variations  in  P-  and  5-wave  speeds  are  resolved  over  an  area  of  approximately  700 
by  400  km  in  central  and  southern  Tibet.  This  result  provides  an  independent  com  parison 
of  the  velocity  structure  beneath  the  Hi-CLIMB  array. 


9.  CONCLUSIONS 


We  have  investigated  frequency-dependent  characteristics  of  Pn  propagation  over  a 
large  region  in  western  China  based  on  nu  merous  events  record  ed  by  the  Hi-CLIMB 
array  from  this  region.  Attributes  of  seismic  wave-trains,  including  arrival  times,  Hilbert 
envelope  amplitudes,  and  instantaneous  and  sp  ectral  frequencies,  have  been  m  odeled  to 
investigate  how  structures  in  the  crust  and  the  upper  mantle  affect  the  propagation  of  the 
Pn  phase  in  Tibet.  For  the  m  odeling,  we  applied  3D  ray  methods,  as  well  as  the  spectral 
element  method  which  provides  complete  numerical  solutions  in  laterally  varying  media. 
The  objective  was  to  obtain  se  lf-consistent  models  of  Pn  propagation  in  western  China 
which  are  free  from  assumptions,  such  as  frequency-independent  geometric  spreading.  To 
this  end,  our  results  should  advance  efforts  in  isolating  effects  of  frequency-dependent 
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propagation  from  those  of  pure-inelasti  c  attenuation  (Q),  1  eading  to  im  proved 
methodologies  for  discrimination  and  yield  estimates  at  regional  distances. 

We  first  modeled  P-wave  arrival-times  to  constrain  the  velocity  stru cture  in  the 
crust  and  the  upper  m  antle  in  central  and  western  Tibet.  T  ravel-times  from  four  events 
along  the  Hi-CLIMB  array  provided  details  on  crustal  velocities,  and  six  events  at 
regional  distances  to  the  array  provide  furt  her  constraints  on  Moho  structure  and  upper 
mantle  lid  velocities.  We  used  3D  ray  trac  ing  to  model  the  travel-times,  and  the  results 
indicated  that  both  Moho  stru  cture  and  m  antle-lid  velocities  in  the  region  are  three- 
dimensional  in  nature  but  approximately  follow  the  trend  of  the  Bangong-Nujiang  Suture 
(BNS).  Although  only  a  limited  number  of  events  were  used  for  the  travel-time  modeling, 
the  results  are  con  sistent  with  earlier  results  from  teleseismic  imaging  using  the  Hi- 
CLIMB  array. 

We  then  modeled  Pn  wave  attributes,  including  e  nvelope  amplitudes  and  pulse 
frequencies,  recorded  by  the  Hi-Climb  array  in  Tibet.  Three  2D  P-wave  velocity  slices 
for  the  events  were  used  in  the  SEM  calculations.  All  the  slices  were  then  earth- flattened 
for  use  in  the  2D  SEM  calcu  lations.  The  SEM  modeling  of  the  pulse  am  plitudes  and 
frequencies  suggested  a  variable  upper  m  antle  structure  in  w  hich  the  Qiangtang  terrane 
has  higher  attenuation  than  the  Lhasa  terrane  in  the  upper  mantle.  However,  the  modeling 
with  a  constant  upper  mantle  attenuation  structure  also  shows  that  this  can  provide  a  good 
fit  to  the  ob  served  amplitudes  and  frequencies,  but  only  if  a  variable  velocity  gradient 
structure  exists  with  a  lower  velocity  gradient  beneath  the  Qiangtang  terrane  compared  to 
the  Lhasa  terrane. 

Thus,  two  alternative  models  for  the  region  are  found,  one  with  a  constant  upper 
mantle  velocity  gradient  similar  to  those  of  Griffin  et  al.  (2011)  and  Phillips  et  al.  (2007) 
with  a  variable  attenuation  m  odel  and  one  with  a  constant  a  ttenuation  model  and  a 
variable  upper  mantle  velocity  gradient  similar  to  that  of  Myers  et  al.  (2010).  In  either 
case,  the  upper  m  antle  of  the  Qiangtang  te  rrane  can  be  distinguished  from  the  Lhasa 
terrane.  The  lower  Pn  velocities  and  lower  Pn  amplitudes  suggest  hotter  temperatures  in 
the  upper  mantle  beneath  the  Qiangtang  terrane  com  pared  to  the  Lhasa  terrane.  Either 
lower  Qp  values  or  a  lower  upper  m  antle  velocity  gradient  in  the  Qiangtang  terrane 
compared  to  the  Lhasa  terran  e  can  be  used  to  match  the  seismic  attribute  data  presented 
here. 

10.  RECOMMENDATIONS 

Since  structural  velocity  m  odels  can  ha  ve  substantial  eff  ects  on  the  frequency 
dependent  propagation  of  Pn  wave  s,  it  is  recommended  that  these  effects  be  further 
studied  using  adequately  sampled  seismic  array  data  from  regional  events.  This  then  can 
be  used  to  establish  baselines  for  specific  locations  for  discrimination  and  yield  estimates 
at  regional  distances.  However,  more  work  is  needed  to  investig  ate  trade-offs  between 
spatially  varying  attenuation  m  odels  and  complex  velocity  m  odels.  Along  with  more 
complicated  structural  models,  frequency  variable  Q  models  from  small-scale  scattering 
also  need  to  be  incorporated  into  the  modeling  of  Pn  waves. 
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APPENDIX  A.  CALCULATION  OF  SEISMIC  ATTRIBUTES 


Let  s(t )  be  a  seismic  wavelet,  then  the  analytical  signal  a(t)  can  be  written  as 


a(t)=s(t)+  is(t ) 


(A.l) 


where 


«o=-J 

TT  J 


1  r  s(j)dr 


7Z  J  T-t 

—00 


(A.2) 


is  the  Hilbert  transform  (Bracewell,  2000).  The  analytical  si  gnal  in  the  time  domain  can 
also  be  written  as 

a(t)=eim  (A.  3) 


where  the  instantaneous  amplitude  (trace  envelope)  is  given  by 


(A.4) 


and  the  instantaneous  phase  is 


(*(0=tair'(^). 

s(t) 


(A.5) 


The  instantaneous  frequency  is  then  given  by 


In  at 


(A.6) 


This  can  also  be  calculated  as 


m= 


1  dt  dt 

In  e 2  (t)  +  s2 


(A.  7) 


where  damping,  e 2 ,  is  used  to  stable  the  instantaneous  frequency  when  the  envelope 
amplitude  is  small.  This  is  less  im  portant  here  since  the  instantaneous  frequencies  are 
evaluated  at  the  peak  of  the  enve  lope  amplitude  of  the  extracted  Pg  and  Pn  pulses.  In 
addition,  as  the  windo  wed  average  of  the  instantaneous  frequency  weighted  by  the 
squared  envelope  becom  es  larger,  the  estim  ate  approaches  the  centro  id  of  the  p  ower 
spectrum,  where  the  frequency  in  tegral  is  over  the  po  sitive  frequencies  (Nowack  and 
Stacy,  2002;  Matheney  and  Nowack,  1995). 

An  example  of  a  data  trace  from  the  Hi-Climb  array  and  a  windowed  trace  an  d 
envelope  are  shown  in  figure  A.l.  For  the  pulse,  the  amplitude  is  taken  at  the  peak  of  the 
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Hilbert  envelope.  The  instantaneous  frequency  is  also  taken  at  the  time  where  the  peak  of 
Hilbert  envelope  occurs.  The  centroid  frequenc  y  is  taken  from  the  Fourier  transform  of 
the  windowed  seismic  pulse.  If  the  envelope  doesn't  suffici  ently  drop  after  the  first 
arriving  pulse,  the  trace  is  not  included  fo  r  further  am  plitude  and  frequency  analysis 
(Matheney  andNowack,  1995).  A  Iso,  ifthees  timated  instantaneous  f  requency  is  not 
consistent  with  the  adjacent  traces  in  distan  ce,  then  the  trace  is  not  included  under  the 
assumption  that  the  pu  lse  has  been  contaminated  with  later  arrivals.  Both  the  observed 
and  the  calculated  data  are  analyzed  using  the  same  processing  steps. 

Estimation  of  Hilbert  Envelope  Amplitude 
And  Instantaneous  Frequency 


Figure  Al.  This  shows  an  example  of  a  data  trace  from  the  Hi-CLIMB  array  with  a  window 
function  at  the  top  and  a  windowed  trace  and  envelope  at  the  bottom. 
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APPENDIX  B.  VALIDATION  OF  THE  SPECTRAL  ELEMENT  METHOD  FOR 
PN  WAVES 

In  order  to  test  the  s  pectral-element  method  (SEM),  com  parisons  are  m  ade 
between  wide-angle  seism  ic  results  obtained  using  SEM  a  nd  the  reflectivity  results  of 
Braile  and  Smith  (1975),  as  well  as,  between  the  SEM  and  the  asy mptotic  results  of 
Cerveny  and  Ravindra  (1971).  For  these  comp  arisons,  in  addition  to  the  synthetic 
seismograms,  various  seismic  attributes  of  the  Pn  are  derived,  including  pulse  amplitudes 
and  pulse  frequencies.  The  SEM  calculations  are  implemented  using  the  parallel  code 
SPECFEM2D  (Komatitsch  and  Vilolette,  1998;  Komatitsch  et  al.,  2005). 

We  first  compare  the  SEM  with  reflectiv  ity  calculations.  For  various  wide-angle 
applications,  the  reflectivity  m  ethod  has  been  used  extensively  to  compute  synthetic 
seismograms  in  vertically  varying  media.  For  exam  pie,  Braile  and  Sm  ith  (1975) 
computed  synthetic  seismogram s  using  the  re flectivity  method,  and  then  examined  the 
effects  of  various  v  elocity-depth  models  on  th  e  calculated  seismic  waves.  One  of  the 
results  of  the  reflectivity  ca  lculations  was  the  determination  of  wide-angle  reflection 
arrivals  from  the  Moho,  PmP,  and  head- wave  arrivals,  Pn. 

The  velocity  gradient  below  the  Moho  is  an  i  mportant  parameter  for  the 
characterization  of  the  Earth's  upper  m  antle  since  it  greatly  affects  Pn  amplitudes.  In 
addition,  the  flattening  of  a  sphe  rical  Earth  for  calculation  purposes  will  resu  It  in  an 
effective  positive  velocity  grad  ient  below  the  Moho.  As  a  resu  It,  an  effective  pos  itive 
velocity  gradient  is  prevalent.  A  p  ositive  velocity  gradient  results  in  diving  wave  s  and 
interference  head-waves,  while  a  negativ  e  velocity  gradient  s  trongly  decays  the 
amplitudes  of  head- waves  (Hill,  1971). 

To  make  a  correct  comparison  between  th  e  SEM  and  the  reflectivity  results,  the 
pulse  shapes  of  the  initial  sources  for  each  m  ethod  must  be  the  same.  In  the  reflectivity 
calculations  of  Braile  and  Smith  (1975),  the  source  time  function  is  given  by 

sin(dt)  —  —sin (mSt)  where  0  <  t  <  t 

„  v  m 

s(t)  =  (B.l) 

0  where  t  <  0  and  t  >  r  , 


with 


5= 


N  7Z 


T 


m=- 


N  +  2 
N 


(B.2) 


where  N  is  an  integer  also  used  to  define  the  number  of  extrem  a  of  the  pulses,  where 
Braile  and  Smith  (1975)  used  N  =  2  .  r  is  the  duration  in  seconds  of  the  source  wavelet, 
and  t  is  the  sampled  time  (Fig.  B.l),  where  Braile  and  Sm  ith  (1975)  used  r=0.2  sec  .  In 
the  SEM  modeling,  the  first  derivative  of  a  Ga  ussian  distribution  is  chosen  to  provide  a 
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similar  source  pulse  sh  ape  (Fig.  B.l).  The  initia  1  source  tim  e  function  intheSEM 
calculations  is  given  by 

s(t)  =  -  2 A(  n2a\  )  e"*2 a°2('“'o)2  (B.3) 

where  A  is  a  multiplication  factor,  a0  is  related  to  the  width  of  the  source  tim  e  function, 
and  t0  is  the  tim  e  delay  (Fig.  B.l).  In  order  to  m  ake  these  pulses  sim  ilar,  a0  =  1 A  is 
found  to  be  appropriate  for  the  SEM  source  time  function. 


Source  Time  Functions  for  Reflectivity  and  SEM 


Figure  B.l.  Source  time  functions  for  the  reflectivity  calculation  (solid  line)  with  r  =  0.2  sec 
and  the  SEM  calculation  with  ao  =7.4  (dashed  line). 

Amplitudes  of  the  seismograms  vary  in  a  different  manner  for  2D  and  3D  seismic 
wave  propagation  problem  s  (Chew,  1990).  Reflectivity  method  amplitudes  from  Braile 
and  Smith  (1975)  are  calculated  for  3D,  wh  ile  SEM  amplitudes  are  calculated  for  2D 
wave  solutions,  and  the  amplitude  variations  with  distance  are  different.  In  the  2D  SEM, 
the  source  can  be  con  sidered  as  a  line  sou  rce.  Since  the  problem  has  cylind  rical 
symmetry,  this  can  b  e  written  in  cylind  rical  coordinates  as  U(x,y)  =  U(R )  and 
S(x) S(y)=S(R)  (Chew,  1990),  where  R  is  the  distance  from  the  source.  Assuming  that 
U (R)  has  a  time  dependence  e~im ,  amplitudes  far  from  the  source  can  be  written  as 


U(R)~C 


ikR 


]j  i  n  kR 


(B.4) 


where  C=H 4,  k=colal,  and  ax  is  the  wave  speed.  Therefore  2D  amplitudes  can  be 
written  as 
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~^—e‘ 

inkR 


(B.5) 


ikR 


By  substituting  i  with  exi2  and  Vl77  with  e  ,ff/4in  Eqn  (B  .5),  we  can  obtain  the 
following  equation  for  2D  amplitude  variations  in  the  frequency  domain  in  the  far  field  as 

.  (B.6) 

An  important  point  here  is  that  the  exponential  term  e,mRla'  in  the  equations  above  results 
in  a  delay  of  the  pulse  in  time.  Because  Matlab  routines  are  used  to  make  all  corrections 
and  Matlab  uses  a  different  sign  convention  in  the  frequency  domain  than  Chew  (1990), 
in  terms  of  a  time  delay  for  a  pulse,  these  exponential  terms  must  have  negative  sign  in 
Matlab  implementation. 

In  the  3D  reflectivity  m  ethod,  the  source  can  be  considered  as  a  point  source. 
Chew  (1990)  defines  a  point  source  in  the  frequency  domain  as 

U(R)  =  —^—  eikR  =—^—  e,taR,a' .  (B.7) 

47rR  4  7rR 

In  order  to  correct  the  2D  SEM  amplitudes  for  3D  reflectivity,  one  needs  to  multiply  the 
2D  SEM  amplitudes  by 


n 


n  co 
2a, R 


(e“!"/4). 


(B.8) 


In  a  Matlab  im  plementation,  this  correction  has  to  be  done  with  a  reverse  sign  on  the 
exponential. 


-JAA  («'"*)  (Matlab).  (B.9) 

7T  \  2  axR 

In  Matlab,  multiplying  the  2D  SEM  amplitudes  by  the  term  in  Eqn.  (B.9),  we  ob  tain  the 
corrected  amplitudes  of  2D  SEM  f  or  comparison  with  the  3D  ref lectivity.  After  this 
correction  is  done,  am  plitudes  of  both  m  ethods  are  multiplied  by  distance  for  plotting 
purposes.  Although,  this  is  correct  for  horn  ogeneous  media,  it  will  still  correct  th  e  2D 
calculations  for  approximate  out  of  plane  sp  reading  for  a  medium  only  varying  in  a  2D 
plane.  Higher  order  asymptotic  corrections  fr  om  2D  to  effective  3D  are  given  by  Miksat 
et  al.  (2008)  as  well  as  other  work  by  Yidale  et.  al  (1985),  Helmberger  and  Vidale  (1988) 
and  Chen  et.  al  (2007). 
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Synthetic  Seismograms  for  Model  Moho-1  (Vertical  Components)  Original  SEM  Amplitudes 
with  no  Corrections  (top)  and  with  Pulse  Shape  and  Corrections  (bottom) 


Figure  B.2.  SEM  synthetics  of  model  Moho-1  with  no  3D  correction  (top)  and  SEM 
synthetics  of  model  Moho-1  with  3D  distance  and  pulse  shape  corrections  (bottom). 
Top  and  bottom  amplitudes  are  scaled  by  the  maximum  amplitude  of  all  traces. 

The  first  com  parison  is  for  a  con  stant  velocity  30  km  crust  over  a  constant 
velocity  mantle,  and  is  called  Moho-1  by  Bra  ile  and  Smith  (1975).  The  P-wave  velocity 
in  the  crust  ax  is  6.4  km/sec  and  the  P-wave  velocity  in  the  mantle  a2  is  8.0  km/sec.  The 

S-wave  velocities  /?  of  the  layers  are  computed  assuming  to  a  Poisson  ratio  =  0.25.  The 

densities  are  computed  from  p;  =  0.252  +  0.3788a;  (Birch,  1964  and  Braile  and  Smith, 
1975).  In  figure  B.2,  the  SEM  calculations  fo  r  Moho-1  are  shown.  In  the  top  plot,  the 
SEM  synthetics  are  scaled  by  the  m  aximum  amplitude  of  all  traces.  The  bottom  plot 
shows  the  SEM  synthetics  multiplied  by  the  Eqn  (B.9)  in  th  e  frequency  domain  for  a  3D 
amplitude  correction,  th  en  scaled  b  y  the  maximum  amplitude  of  all  traces.  Thes  e  two 
plots  show  how  the  pulse  shapes  vary  between  2D  and  3D  calculations.  In  figure  B.3,  a 
plot  of  the  results  of  Braile  and  Sm  ith  (1975)  for  Moho-1  (solid  line)  and  those  from  the 
corrected  results  of  SEM  (dashed  line)  are  shown.  Both  results  are  multiplied  by  distance 
for  plotting  purposes.  As  can  be  seen  from  the  comparison,  the  resu  Its  are  in  good 
agreement. 
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Synthetic  Seismograms  for  Model  Moho- 1 


Figure  B.3.  A  comparison  between  the  SEM  calculations  (dashed)  and  the  reflectivity 
results  (solid)  of  Braile  and  Smith  (1975)  for  a  constant  velocity  30  km  crust  over  a 
constant  velocity  mantle.  The  amplitudes  of  the  reflectivity  results  and  the  3D 
corrected  SEM  results  are  both  multiplied  by  distance  in  the  plot. 

The  second  comparison  model  is  a  consta  nt  velocity  30  km  crust  over  a  positiv  e 
velocity  gradient  in  the  top  5  km  of  the  mantle,  and  is  nam  ed  Moho-5  by  Braile  and 
Smith  (1975).  The  P-wave  velocity  of  the  crust  ax  is  6.4  km  /sec  which  is  the  sam  e  as 
Moho-1.  However,  the  velocity  gradient  in  the  upper  m  antle  for  5  km  below  the  Moho  is 
calculated  by  a2(z )  =  a^p  +  bz  where  a^p  is  the  top  velocity  of  the  mantle  (8  km/sec), 
b  is  the  velo  city  gradient  (0.08  sec-1),  and  z  is  the  depth.  At  the  depth  of  35  km,  the 
velocity  of  the  mantle  is  8.4  km /sec  and  this  velocity  is  constant  down  to  the  bo  ttom  of 
the  model  (60  km).  The  computation  of  the  S-  wave  velocities  and  de  nsities  are  done  in 
the  same  manner  as  the  Moho-1  model  shown  above.  Similar  to  the  Moho-1  m  odel,  the 
comparisons  between  the  reflectivity  results  of  Braile  and  Smith  (1975)  and  those  from 
the  corrected  SEM  results  are  in  good  agreement  (figure  B.4). 

For  the  Moho-1  model,  figure  B.3  shows  that  sub-critical  reflections  (less  than  80 
km  in  distance)  reproduce  the  initial  s  ource  time  function.  For  the  region  around  the 
critical  distance,  because  t  he  reflected  waves  and  head- waves  begin  interfering,  larger 
amplitude  arrivals  start  appearing.  There  is  a  zone  around  the  critical  distance  that  larger 
amplitudes  can  be  observed.  Th  e  size  of  th  is  zone  depends  on  the  frequency  content  of 
the  source  and  lower  frequencies  will  cause  a  larger  sized  zone  (Braile  and  Smith,  1975). 
In  figure  B.3  and  B.4,  wide-angle  reflections  can  easily  be  observed  sin  ce  they  have  the 
largest  amplitudes  compared  to  th  e  other  arrivals.  Recalling  that  the  travel  tim  es  are 
reduced  8  km/sec  and  by  looking  at  the  given  velo  city  slope  chart  in  the  figures,  it  can  be 
seen  that  wide-angle  reflections  have  an  apparent  velocity  of  6.4  km  /sec  which  is  the 
upper  homogenous  layer's  velocity.  In  figures  B.3  and  B.4,  the  Pn  amplitudes  for  the 
Moho-5  model  are  much  larger  than  the  head- wave  amplitudes  for  the  Moho-1  model. 
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Synthetic  Seismograms  for  Model  Moho-5 
(Vertical  Component) 


Figure  B.4.  A  comparison  between  the  results  from  the  SEM  calculations  (dashed)  and  the 
reflectivity  results  (solid)  of  Braile  and  Smith  (1975)  for  a  constant  velocity  30  km 
crust  over  a  mantle  with  a  velocity  gradient  in  the  top  5  km.  Amplitudes  of  the 
reflectivity  results  and  the  corrected  SEM  results  are  both  multiplied  by  distance  in 
the  plot. 

This  effect  is  due  to  the  presence  of  a  positive  velocity  gradient  below  the  interface  (Hill, 
1971).  In  figure  B.4,  interference  head-waves  start  appearing  at  distances  where  the  pure 
head- waves  are  observed.  In  general,  interference  head- waves  have  very  sim  ilar  travel¬ 
time  characteristics  to  the  pure  head-waves,  but  curve  for  greater  distances.  However,  the 
amplitude  characteristics  of  interference  head-waves  are  q  uite  different  from  the  pure 
head- wave  amplitudes  (Braile  and  Smith,  1975). 

SEM  is  next  com  pared  with  the  asym  ptotic  results  of  Cerveny  and  Ravindra 
(1971).  Only  one  interface  is  considered,  and  source  and  th  e  receivers  are  located  in  the 
first  medium.  The  pure  head  wav  e  case  of  Ce  rveny  and  Ravindra  is  equivalent  to  the 
Moho-1  model  given  above  of  Braile  and  Sm  ith  (1975).  The  pure  head  wave  case  has  a 
constant  velocity  30  km  crust  over  a  consta  nt  velocity  m  antle  model.  Also,  P-wave 
velocity  in  the  crust  is  6.4  km/sec  and  in  the  mantle  is  8  km/sec  which  are  identical  to  the 
Moho-1  model  of  Braile  and  Sm  ith  (1975).  Cerveny  and  Ravindra  (1971)  chose  the 
densities  of  different  layers  to  be  the  same  pxl  p2=\.  The  interference  head  wave  case  of 
Cerveny  and  Ravindra  (1971)  is  a  constan  t  velocity  30  km  crust  over  a  positiv  e  velocity 
gradient  mantle.  The  velocity  in  the  lower  medium  is  given  by  a2(z)  =  a2°P  +  bz  to 
calculate  the  velocity  gradient  where  a 2°p  is  the  top  velocity  of  the  mantle  (8  km/sec)  b  is 
the  velocity  gradient  (0.002264  sec'1),  and  z  is  the  depth. 
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Amplitude  vs  Distance 


Distance  (km) 


Figure  B.5.  This  shows  a  comparison  between  the  SEM  calculations  (crosses)  and  the 
asymptotic  calculations  (background  plot)  of  Cerveny  and  Ravindra  (1971)  for  the 
pure  and  interference  head  waves.  C+  wave  is  the  interference  head  wave,  C0  is  the 
diving  wave,  and  the  131  wave  is  the  pure  head  wave  case  (Modified  background 
plot  from  Cerveny  and  Ravindra,  1971). 


There  is  good  agreement  between  asym  ptotic  calculations  of  Cerveny  and 
Ravindra  (1971)  and  the  SEM  results  for  pure  and  interference  head- wave  cases.  A 
positive  velocity  gradient  below  the  in  terface  drastically  increases  the  h  ead-wave 
amplitudes  (Cerveny  and  Ravindra,  1971).  F  igure  B.5  shows  the  comparison  of  the 
amplitudes  of  pure  and  interferen  ce  head  wa  ves  for  the  asym  ptotic  calculations  of 
Cerveny  and  Ravindra  (1971)  (background  plot)  and  the  SEM  results  (crosses).  The  pure 
head-wave  case  of  Cerveny  and  Ravindra  (1971)  is  referred  to  by  131  in  figure  B.5.  Also 
the  interference  head-wave  case  of  Cerveny  and  Ravindr  a  (1971)  is  referred  to  by  C  +  in 
figure  B.5.  Because  the  amplitude  variations  with  distan  ce  are  d  ifferent  for  the  3D 
asymptotic  method  and  2D  SEM  (Chew,  1990)  as  mentioned  above,  the  same  amplitude 
correction  in  frequency  domain  has  been  made  here  for  the  2D  SEM  results  to  com  pare 
with  the  3D  asymptotic  results  of  Cerveny  and  Ravindra  (1971). 

The  amplitude  spectrum  as  a  function  of  distance  is  given  in  figure  B.6.  The 
initial  pulse  has  a  dom  inant  frequency  of  6.4  Hz.  The  initia  1  lowering  of  the  ce  ntroid 
frequency  of  the  spectra  is  at  140  km.  Then  the  centroid  frequency  starts  increasing.  It 
then  gets  higher  than  the  incident  wave  centroid  frequency  at  180  km  .  For  greater 
distances,  the  diving  waves  begin  to  dom  inate.  The  em  ergence  of  the  diving  waves 
causes  a  lowering  of  the  centroid  frequency  m  aking  it  s  imilar  to  the  incident  wave 
centroid  frequency  again.  This  is  consiste  nt  with  the  asym  ptotic  Gaussian  b  earn 
calculations  of  Nowack  and  Stacy  (2002).  Howe  ver,  there  is  a  tuning  effect  between 
different  interference  wave  types,  and  th  is  can  result  in  an  overshoot  to  higher 
frequencies  for  the  centroid  frequency. 
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Interference  Head  Wave  Model 


Figure  B.6.  This  shows  amplitude  spectra  as  a  function  of  distance  for  the  interference  head 
wave  case.  The  thin  black  curves  are  from  Cerveny  and  Ravindra  (1971)  and  the 
thicker  curves  are  from  SEM  calculations  for  the  same  model  (Background  plot 
from  Cerveny  and  Ravindra,  1971). 

The  SEM  is  finally  compared  with  an  analytical  solu  tion  for  an  attenuating 
medium.  Amplitudes  and  pulse  shapes  in  an  attenuating  medium  are  different  from  the 
amplitudes  in  an  elastic  m  edium,  since  amplitudes  decay  expone  ntially  due  to  the 
attenuation.  Aki  and  R  ichards  (2002)  give  th  e  amplitude  decay  in  frequency  fo  r  an 
attenuating  medium  with  a  constant  Q  as; 

coR 

U(x,a>)  =  U0(co)  e  2aQ  (B.10) 


where  co=2n  f ,  R  is  the  dis  tance,  a  is  the  velocity  of  the  medium,  Q  is  the  qu  ality 
factor  of  the  medium,  and  U(x,a> )  is  the  attenuated  amplitude.  Note  that  this  attenuation 
filter  will  give  an  acausal  operator  and  an  added  phase  term  is  needed  for  a  causal  filter 
(Aki  and  Richards,  2002).  However  this  is  su  itable  for  testing  amplitudes  and  amplitude 
spectral  shapes.  To  test  the  attenua  tion  of  the  SEM,  a  homogeneous  viscoelastic  model 
with  a  Qp  of  1 12.5,  a  P-wave  veloc  ity  of  5  km/sec,  and  a  dominant  frequency  of  2  Hz  is 

calculated  for  various  s  ource  receiver  distances.  The  synth  etic  traces  of  the  SEM  are 
compared  with  the  analytical  calculations  for  the  same  Q  value  applied  in  the  Eqn  (B.10). 
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SEM  Waveforms  Compared  to  Attenuated  Waveforms 

Qp  =  1 12.5 


Source  at  1 0  km 
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Figure  B.7.  This  shows  the  comparison  of  the  SEM  traces  with  the  analytical  calculations 
for  a  homogeneous  viscoelastic  model  with  a  Qp  of  112.5. 


For  this  test,  a  receiver  is  put  at  th  e  surface  and  four  shots  at  different  depths 
vertically  below  the  receiver  are  us  ed  to  calculate  the  SEM  synthetics.  Then,  the  SEM 
synthetics  are  compared  with  a  con  stant  Q  (acausal)  filter  applied  to  the  near  distance 
SEM  trace  at  10  km  depth.  For  this  test,  the  bulk  attenuation  is  neglected  by  using  a  large 
Qk  value  (  Qk  =1000 )  and  the  Poisson's  ratio  is  assumed  to  be  0.25  (  a2  /  /?2  =  3 )  in  Eqn 
(B.l  1)  (from  Dahlen  and  Tromp,  1998). 


Q, 


=(i- 


4  f52 
3  a1 


a 


A4)- 

3  V  a 


(B.l  1) 


where  a  and  [1  are  the  P  an  d  S  wave  velocities  of  the  m  edia.  Then  Q ,  equal  to  Qs ,  is 
chosen  to  be  50  resulting  in  a  Qp of  1 12.5  from  the  ratio  of  Qp  =  9  Qs  /  4  between  Qp  and 
Qs  for  a  larg  e  QK.  In  orde  r  to  reta  in  this  ratio,  all  S  EM  calculations  in  the  f  olio  wing 

chapters  are  made  based  on  assum  ing  attenuation  loss  in  shear.  Figure  B.7  shows  the 
original  SEM  traces  compared  with  the  analytically  calculated  traces.  Because  an  acausal 
Q  model  was  used,  the  analytic  traces  were  aligned  with  the  SEM  traces  to  highligh  t  the 
effects  of  pulse  broad  ening  and  am  plitude  decay.  Figure  B.8  shows  the  envelope 
amplitudes  of  the  original  SEM  outputs  com  pared  with  the  envelope  amplitudes  of  the 
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analytically  calculated  traces.  Figure  B.9  shows  the  frequency  spectrum  comparisons  of 
the  original  SEM  traces  with  the  analytically  calculated  traces. 


Amplitudes  vs  Distance  Qp  =  1 12.5 
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Figure  B.8.  This  shows  the  envelope  amplitude  of  the  SEM  traces  with  the  envelope 
amplitude  of  the  analytically  calculated  traces  given  in  figure  B.7. 
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Figure  B.9.  This  shows  a  comparison  of  the  frequency  spectra  of  the  SEM  traces  with 
analytical  calculations  for  the  traces  given  in  figure  B.7. 


46 


In  this  appendix,  synthetic  seismograms  using  the  SE  M  approach  are  compared 
with  the  reflectivity  results  from  Braile  and  Smith  (1975)  and  the  asymptotic  calculations 
of  Cerveny  and  Ravindra  (1971).  T  hese  comparisons  show  that  th  e  SEM  approach  can 
handle  homogeneous,  inhomogeneous  and  viscoelastic  models,  where  very  similar  results 
are  obtained  from  the  SEM  compared  to  the  other  approaches.  However  a  correction  for 
out-of-plane  spreading  is  needed  to  com  pare  2D  and  3D  results  in  the  case  of  only  2D 
variations  in  the  medium  as  used  here.  For  a  positive  velocity  gradient  below  the  Moho, 
an  interference  head- wave  will  result,  and  the  amplitude  and  frequency  characteristics 
can  be  quite  different  than  the  pure  head-wave  case. 
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Abstract 


Using  data  from  regional  earthquakes  recorded  by  the  Hi-CLIMB  array  in  Tibet, 
we  model  P-wave  arrival-times  to  constrain  the  velocity  structure  in  the  crust  and  the 
upper  mantle  in  central  and  western  Tibet.  Of  more  than  30  high-quality,  regional  seismic 
profiles  that  have  been  assembled,  we  have  selected  10  which  show  excellent  crustal  and 
Pn  arrivals  for  further  analysis.  Travel  times  from  four  events  along  the  Hi-CLIMB  array 
provide  details  on  crustal  velocities,  and  six  events  at  regional  distances  to  the  array 
provide  further  constraints  on  Moho  structure  and  upper  mantle  lid  velocities.  We  use  3- 
D  ray  tracing  to  model  the  travel  times,  and  the  results  indicate  that  the  Moho  beneath  the 
Lhasa  terrane  of  southern  Tibet  is  over  73  km  deep  with  a  Pn  speed  of  about  8.2  km/s. 
The  Qiangtang  terrane  north  of  the  Bangong-Nujiang  suture  (BNS)  shows  a  thinner  crust, 
by  up  to  10  km,  and  a  lower  Pn  speed  of  7. 8-7. 9  km/s.  Travel  times  from  events  to  the 
west  and  east  of  the  array  indicate  that  both  Moho  structure  and  mantle-lid  velocities  in 
the  region  are  three-dimensional  in  nature  but  approximately  follow  the  trend  of  the  BNS. 
Although  only  a  limited  number  of  events  were  used  for  the  travel  time  modeling,  the 
results  are  consistent  with  earlier  results  from  teleseismic  imaging  using  the  Hi-CLIMB 
array. 

Introduction 

The  Himalayan-Tibetan  orogeny  is  the  largest  active  continent-continent  collision 
zone  in  the  world  (Yin  and  Harrison,  2000),  and  the  continuing  convergence  has  resulted 
in  large  scale  deformation  in  central  Eurasia  (Molnar  and  Tapponnier,  1975;  Tapponnier 
et  al.,  2001).  The  Hi-CLIMB  (Himalayan-Tibetan  Continental  Lithosphere  During 
Mountain  Building)  experiment  was  conducted  along  a  north-south  corridor  between 
84°E  and  86°E  extending  from  the  Ganga  foreland  basin,  over  the  Lesser  and  Higher 
Himalayas,  crossing  the  Indus- Yarlung  suture  and  the  Bangong-Nujiang  suture,  then 
reaching  into  the  Qiangtang  terrane  in  central  Tibet  (Figure  1).  Between  2002  to  2005, 
the  Hi-CLIMB  array  was  deployed  in  3  phases  with  a  total  of  over  210  broadband  sites 
(Nabelek  et  al.,  2005).  A  unique  feature  of  the  Hi-CLIMB  experiment  is  that  it  crosses 
the  highest  portions  of  the  Himalayas  from  the  Indian  shield  to  the  interior  of  Tibet.  It 
also  has  a  unique  combination  of  a  dense  station  spacing  of  3  to  8  km  for  the  main  linear 
array  and  a  broad  regional  coverage.  In  this  study  we  investigate  data  from  the  linear 
array  component  of  Hi-CLIMB  in  Tibet  to  analyze  first  arrival  travel  times  from  regional 
earthquakes. 

Because  of  the  importance  of  the  region  for  understanding  active  continental 
collision,  Tibet  has  been  the  site  of  several  large  scale  field  experiments  in  recent  years 
(e.g.  Him  et  al.,  1984;  Zhao  et  al.,  1993;  McNamara  et  al.,  1994;  Kind  et  al.,  2002; 
Tilmann  et  al.,  2003;  Wittlinger  et  al.,  2004,  2009).  Studies  of  Pn  in  Tibet  have  found 
thicknesses  of  roughly  double  the  continental  average  (Chen  and  Molnar,  1981;  Heam  et 
al.,  2004;  and  Liang  et  al.,  2004,  2006;  Sun  and  Toksoz,  2006). 

Using  teleseismic  data  from  the  Hi-CLIMB  experiment,  the  crust  was  found  to  be 
between  70  and  80  km  thick  beneath  the  Lhasa  terrane,  but  on  the  order  of  10  km  thinner 
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beneath  the  Qiangtang  terrane  to  the  north  (Nabelek  et  ah,  2009;  Tseng  et  al,  2009; 
Nowack  et  ah,  2010).  From  the  analysis  of  wide-angle  and  teleseismic  data  from  the 
INDEPTH  III  experiment  located  about  5  degrees  to  the  east  of  Hi-CLIMB  but  crossing 
similar  terranes,  a  somewhat  smaller  change  in  crustal  thickness  was  found  from  south  to 
north  (Zhao  et  ah,  2001;  Kind  et  ah,  2002;  Shi  et  ah,  2004,  Tian,  et  ah,  2005).  While  the 
Moho  conversion  is  distinct  and  coherent  beneath  a  majority  of  both  the  Lhasa  and 
Qiangtang  terranes  from  receiver  function  images  of  Hi-CLIMB  data  images,  the 
conversion  strength  was  observed  to  vary  with  the  direction  of  illumination  (Nabelek  et 
ah,  2009;  Nowack  et  ah,  2010).  Also  from  Hi-CLIMB  data,  a  zone  of  crustal  and  Moho 
disruption  marks  the  Bangong-Nujiang  suture  (BNS)  zone  between  the  more  stable  Lhasa 
and  Qiangtang  blocks  (Nowack  et  ah,  2010). 

From  travel  time  studies,  the  average  Pg  velocities  in  the  Tibetan  crust  was  found 
to  be  lower  than  the  continental  average,  in  the  range  of  6.2  to  6.3  km/s  (Zhao  et  ah, 

2001;  Monsalve  et  ah,  2008;  Steck  et  ah,  2009).  Studies  from  the  Hi-CLIMB  and 
INDEPTH  experiments,  which  have  more  dense  station  spacings,  found  similar  average 
crustal  velocities.  Also,  lower  crustal  velocities  were  found  in  the  areas  of  the  Indus- 
Yarlung  suture  (IYS)  and  the  Bangong-Nujiang  suture  (BNS)  for  Hi-CLIMB  (Hung  et 
ah,  2010)  and  in  the  area  of  the  BNS  for  INDEPTH  III  (Haines  et  ah,  2003;  Meissner  et 
ah,  2004). 

Upper  mantle  P- wave  velocities  beneath  Qiangtang  terrane  north  of  the  BNS  were 
found  from  regional  travel  time  studies  to  be  as  low  as  7.9  km/s  (Liang  et  ah,  2004; 

Hearn  et  ah,  2004;  Pei  et  ah,  2007).  These  low  velocities,  as  well  as  the  inefficient  Sn 
propagation  observed  in  this  area  (Barazangi  and  Ni,  1982;  Ni  and  Barazangi,  1983; 
Brandon  and  Romanowicz,  1986;  McNamara  and  Owen,  1995;  McNamara  et  ah,  1997; 
Xie  et  ah,  2004;  Barron  and  Priestley,  2009),  imply  that  the  thinner  crust  of  the 
Qiangtang  terrane  is  at  least  in  part  isostatically  supported  by  a  mantle  beneath  that  is 
hotter  and  more  buoyant  (Zhao  et  ah,  2001;  Jimenez- Munt  et  ah,  2008).  The  P-wave 
velocities  inferred  beneath  the  southern  Tibetan  Plateau,  south  of  the  Bangong-Nujiang 
suture,  are  as  high  as  8.3  km/s  and  similar  to  the  upper  mantle  velocities  observed 
beneath  stable  cratonic  regions,  such  as  the  Indian  shield  or  the  Tarim  and  Sichuan  basins 
(McNamara  et  ah,  1997;  Liang  and  Song,  2006;  Pei  et  ah,  2007). 

The  differences  between  results  from  Hi-CLIMB  and  INDEPTH  III  highlight  the 
fact  that,  while  the  Tibetan  Plateau  is  composed  of  a  succession  of  terranes  from  north  to 
south  (Dewey  et  ah,  1988),  there  exists  significant  lateral  variation  from  west  to  east,  and 
the  structure  of  the  plateau  must  be  taken  into  account  in  all  three  dimensions.  In  this 
study,  forward  modeling  of  P-wave  travel  times  from  regional  earthquakes  recorded  by 
the  linear  array  component  of  Hi-CLIMB  in  Tibet  is  used  to  determine  the  three- 
dimensional  P-wave  velocity  structure  for  the  region  of  the  central  and  western  Tibetan 
Plateau  beneath  the  Hi-CLIMB  array. 
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Data  Analysis 

We  used  broadband  seismic  data  from  regional  earthquakes  in  Tibet  recorded  by 
the  densely  spaced,  linear  component  of  the  Hi-CLIMB  array  in  Tibet  (Figure  1)  to 
investigate  crustal  and  upper  mantle  structure.  Based  on  locations  of  epicenters  reported 
in  the  EHB  catalog  (Engdahl  et  ah,  1998,  see  Data  and  Resources)  or  the  NEIC  PDE 
catalog  (USGS,  See  Data  and  Resources;  for  events  not  included  in  the  EHB),  we 
selected  events  for  which  both  Pg  and  Pn  arrivals  were  recorded  by  the  Hi-CLIMB  array 
in  Tibet  (Figure  1;  a  list  of  events  is  available  as  an  electronic  supplement  to  this  paper). 

To  construct  a  velocity  model  along  the  trend  of  the  Hi-CLIMB  array  in  Tibet, 
we  relied  on  data  from  four  events  located  along  and  nearby  the  array  (L2a,  L2d,  L3a, 
and  L3b;  Figure  1)  for  crustal  structure.  In  addition,  we  used  six  events  whose  epicenters 
are  about  200  to  300  km  away  from  the  nearest  station  to  constrain  seismic  velocity 
structures  in  the  crust  and  the  upper  mantle  on  a  regional  scale.  Four  of  the  events  to  the 
south  of  the  array  have  similar  epicenters.  This  provides  some  verification  on  the 
variation  in  the  P-wave  travel  times  due  to  small-scale  heterogeneity  as  well  as  picking 
errors.  The  entire  group  of  events  was  selected  to  be  largely  in-line  with  the  trend  of  the 
array,  but  also  to  provide  a  range  of  azimuthal  coverage.  We  visually  picked  the  timing 
of  the  first  P-arrivals  on  the  vertical  component  of  the  broadband  seismograms  after 
applying  a  six-pole,  Butterworth  band-pass  filter  between  0.5  and  5.0  Hz.  The  average 
pick  error  is  estimated  to  be  +/-  0.3  s.  Figure  2  shows  examples  of  seismic  record- 
sections  from  two  events  with  the  traces  plotted  with  respect  to  surface  distance  of  the 
event  to  each  station,  and  all  travel  times  are  reduced  by  8.0  km/s. 

We  carried  out  forward  modeling  of  the  travel  times  using  a  3D  ray  tracing  code 
by  Cerveny  et  al.  (1988;  the  CRT  algorithm).  In  doing  so,  the  elevation  of  the  free 
surface  in  the  model  is  at  5  km  above  sea-level,  the  average  elevation  of  the  Hi-CLIMB 
stations.  The  model  consists  of  a  grid  system,  with  a  spacing  of  1°  in  longitude,  0.5°  in 
latitude,  and  10  km  in  depth,  over  the  region  between  27  to  36°N  in  latitude  and  80  to 
90°E  in  longitude  and  extending  down  to  a  depth  of  200  km.  The  Moho  is  an  explicit 
surface  separating  grid-points  in  the  crust  from  those  in  the  upper  mantle.  A  scheme  of 
cubic  splines  interpolate  among  the  nodes  of  the  grid  system  to  specify  values  of  seismic 
wave  speeds  and  positions  of  the  Moho  and  to  calculate  ray  paths  in  the  three- 
dimensional  model.  An  example  of  the  3D  ray  tracing  for  event  Lie  is  shown  in  Figure 
3. 

Earthquake  Events  along  the  Hi-CLIMB  Array 

The  epicenters  for  earthquakes  L2a,  L2d,  L3a,  and  L3b  are  located  along  the 
length  of  the  Hi-CLIMB  array  in  Tibet  (Figure  1),  and  the  Pg  travel  times  recorded  from 
these  nearby  sources  are  used  to  constrain  the  crustal  velocity  structure  beneath  the  array. 
The  four  events  were  first  relocated  using  the  location  code  HYPOINVERSE  2000 
(Klein,  2002),  and  the  relocated  epicenters  range  from  4  to  15  km  distant  from  the 
reported  epicenters  (see  the  electronic  supplement  for  a  table  of  the  relocated  hypocenters 
for  the  nearby  events).  Also,  while  the  reported  depths  of  these  sources  are  between  15 
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and  19  km  below  the  surface,  their  depths  changed  to  6-15  km  below  the  surface  after 
relocation,  which  are  more  consistent  with  the  shallow  crustal  seismogenic  zone  in  Tibet 
(Chen  and  Molnar,  1983;  De  la  Torre  et  al.,  2007). 

For  the  sources  along  the  array,  it  can  be  observed  in  Figure  4  that  the  slopes  of 
the  Pg  branches  for  all  four  events  are  generally  linear  both  to  the  north  and  to  the  south 
of  each  source.  The  apparent  velocities  implied  by  the  slopes  of  the  Pg  branches  are 
between  5.9  and  6.1  km/s.  Modeling  of  the  Pg  travel  times  for  these  events  was  first 
done  using  an  initial  ID  crustal  model  (Figure  4 A,  solid  line)  derived  from  an  average 
crustal  model  of  Meissner  et  al.  (2004)  for  the  INDEPTH  III  experiment,  but  with 
somewhat  slower  velocities  in  the  mid  to  lower  crust  as  inferred  by  the  linear  slopes  of 
the  observed  Pg  arrivals  for  Hi-CLIMB.  Fig  4B  shows  the  observed  P- wave  travel  times 
for  the  relocated  events  L2a,  L2d,  L3a,  and  L3b  compared  with  the  travel  times 
calculated  using  the  initial  ID  crustal  velocity  model.  All  travel  times  are  reduced  by  8.0 
km/s  and  plotted  with  respect  to  the  minimum  recorded  travel  time  rather  than  an 
absolute  origin  time  whose  precise  value  is  unknown.  In  Figs.  4B  and  4C,  the  origin 
position  of  each  event  on  the  x-axis  is  determined  by  its  projected  position  along  the 
linear  trend  of  the  Hi-CLIMB  array.  But,  as  with  the  record  sections,  the  P- wave  travel 
times  are  plotted  in  either  direction  away  from  the  source  position  as  a  function  of 
epicentral  distance  from  the  source  to  each  receiver. 

Although  the  initial  ID  crustal  velocity  model  does  a  reasonable  job  of  fitting  the 
travel  times,  there  are  still  some  mismatches  between  the  observed  and  calculated  travel 
times.  The  misfit  between  the  observed  and  calculated  travel  times  can  be  seen  in  Figure 
4B  for  event  L2a  at  distances  less  than  -400  km  and  greater  than  -100  km  on  the  x-axis, 
for  L2d  at  distances  less  than  -300  km  and  greater  than  +50  km  on  the  x-axis,  and  for  L3b 
at  distances  less  than  -50  km  on  the  x-axis.  This  indicates  that  there  is  some  additional 
lateral  variability  in  the  crustal  velocities  along  the  profile  of  the  array.  Therefore, 
starting  from  the  initial  ID  velocity  profile,  a  refined  crustal  velocity  model  was  obtained 
that  includes  a  reduction  of  up  to  0.15  km/s  or  by  around  3  percent  from  the  initial  ID 
crustal  model  in  the  top  30  km  of  the  crust  to  the  north  of  latitude  33.5°N  and  to  the  south 
of  latitude  31°N. 

Figure  4C  compares  the  observed  and  calculated  travel  times  for  the  refined 
crustal  velocity  model  and  shows  that  the  reduced  velocities  beneath  the  southern  and 
northern  portions  of  the  array  provide  a  better  fit  to  the  observed  travel  times  at  far 
offsets.  The  observed  Pg  travel  times  for  event  L3a  are  well  fit  by  the  calculated  travel 
times  using  the  refined  crustal  velocity  model,  and  the  value  of  RMS  misfit  for  the  Pg 
travel  times  is  only  0.25  s,  smaller  than  average  picking  error.  The  refined  crustal 
velocity  model  also  produces  a  good  fit  between  the  observed  and  calculated  Pg  travel 
times  for  event  L3b,  and  the  average  misfit  of  Pg  arrivals  is  0.29  s.  For  event  L2a,  the 
observed  Pg  travel  times  at  far  offsets  are  well  fit  by  those  calculated  for  the  refined 
crustal  velocity  model  and  results  in  an  RMS  misfit  of  0.24  s  for  crustal  arrivals.  The 
observed  travel  times  to  the  north  of  event  L2d  are  generally  well  matched  by  the 
calculated  travel  times  for  the  refined  crustal  velocity  model,  but  the  observed  Pg  travel 
times  for  stations  to  the  south  are  still  somewhat  later  than  those  calculated.  Also,  the 
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observed  travel  times  of  L2d  are  earlier  than  those  calculated  for  this  model  between 
offsets  of  75  to  200  km  to  the  north  which  could  possibly  be  due  to  the  presence  of  a 
localized  zone  of  high  velocity.  However  this  is  not  observed  in  the  travel  times  for  other 
sources  recorded  at  the  same  stations,  so  while  the  average  value  of  RMS  misfit  of  0.43  s 
is  slightly  larger  than  the  average  picking  error,  a  crustal  model  which  provides  a  better 
fit  for  travel  times  at  these  offsets  for  event  L2d  would  not  be  consistent  with  the  data 
observed  from  the  other  events.  Nonetheless,  the  refined  crustal  velocity  model  provides 
the  best  average  fit  to  the  Pg  travel  times  for  all  the  events  along  the  array. 

While  the  Pn  branches  recorded  for  the  nearby  events  shown  in  Figure  4  are 
typically  short  and  have  a  lower  signal-to-noise  ratio,  the  Pg  to  Pn  crossover  distances 
can  offer  some  initial  constraints  on  the  depth  of  the  Moho  along  the  Hi-CLIMB  array. 
The  Pg  to  Pn  crossover  distances  for  the  events  L2a,  L3a,  and  L3b  are  observed  at 
epicenter  distances  in  the  range  of  250  to  300  km  and  are  consistent  with  crustal 
thicknesses  of  up  to  73  km  in  the  south  beneath  event  L3a  but  with  a  thinner  crust  to  the 
north  beneath  event  L3b.  The  Moho  depths  of  the  refined  crustal  velocity  model  can  be 
constrained  by  the  observed  Pg  to  Pn  crossover  distances  of  these  two  events,  as  well  as 
the  observed  Pn  branches  of  the  more  regional  events  located  farther  from  the  array. 
While  not  all  the  observed  Pn  travel  times  for  event  L2a  are  well  fit,  there  is  some 
uncertainty  in  the  Pn  picks  for  this  event  and  there  is  a  possibility  that  these  arrivals  were 
picked  late  or  on  a  different  phase.  However,  this  could  also  result  from  additional  three- 
dimensional  complications  in  the  derived  Moho  model. 

Earthquakes  at  Far  Regional  Distances  from  the  Hi-CLIMB  Array 

Six  far  regional  earthquakes  recorded  by  the  Hi-CLIMB  array  in  Tibet  were 
selected  to  study  Pn  velocities  and  Moho  depths  to  determine  the  velocity  structure  of  the 
crust  and  upper  mantle  in  the  region.  The  events  shown  on  Figure  1  were  selected  which 
have  prominent  Pn  arrivals,  as  well  as  some  Pg  arrivals  so  that  the  crossover  distance  is 
included.  While  these  events  are  located  too  far  from  the  Hi-CLIMB  array  to  be 
effectively  relocated,  the  depths  were  estimated  to  be  consistent  with  the  observed  Pg  to 
Pn  cross-over  distances  for  each  event,  and  the  Moho  depths  inferred  from  the  nearby 
events  and  other  far  regional  events,  (available  as  an  electronic  supplement  to  this  paper; 
see  also,  Griffin,  2010). 

The  observed  P- wave  travel  times  for  the  regional  earthquakes  are  shown  in 
Figure  5.  In  Figure  5A  the  observed  travel  times  for  event  Lie  show  that  the  Pg  travel 
times  are  consistent  with  the  average  velocities  of  around  6.0  km/s  inferred  from  the  four 
events  nearby  the  Hi-CLIMB  array.  It  can  be  seen  that  the  Pn  travel  times  for  stations  at 
distances  between  250  and  450  km  are  delayed  in  comparison  with  those  at  distances 
greater  than  500  km  to  the  north  of  the  event.  Travel  time  modeling  indicates  that  the 
Moho  is  deeper  in  the  southern  portion  of  the  model  at  a  depth  of  around  73  km  and 
shallower  to  the  north  at  a  depth  of  around  64  km.  The  slope  of  the  Pn  branch  to  the 
south  can  be  modeled  by  upper  mantle  velocities  and  velocity  gradients  of  around  8.3 
km/s  and  0.004  s'1,  respectively.  The  slope  of  the  Pn  branch  to  the  north  can  be  modeled 
by  upper  mantle  velocities  and  velocity  gradients  of  7.8  km/s  and  0.003  s'1,  respectively. 
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The  Pn  travel  time  curve  for  event  Lie  shows  a  transition  between  400  and  550  km  from 
the  source  and  indicates  that  the  transition  between  the  deeper  Moho  in  the  south  and  the 
shallower  Moho  to  the  north  occurs  over  a  horizontal  distance  of  roughly  150  km 
between  the  latitudes  of  31°  and  32.5°N.  The  value  of  the  RMS  misfit  between  the 
observed  and  calculated  travel  times  for  this  event  using  the  inferred  3D  regional  velocity 
model  (Figures  6  and  7)  is  0.27  s.  Some  of  the  misfit  might  suggest  further  complications 
on  the  Pn  speeds  away  from  the  Hi-CLIMB  array.  However,  here  a  self-consistent 
velocity  model  among  all  the  seismic  gathers  is  developed. 

Figure  5B  shows  the  observed  travel  times  for  event  Lla  located  to  the  south  of 
the  Hi-CLIMB  array.  As  with  event  Lie,  the  observed  Pg  travel  times  are  consistent 
with  the  crustal  velocities  observed  using  the  four  nearby  events  to  the  Hi-CLIMB  array. 
The  observed  Pg  to  Pn  crossover  distance  and  Pn  travel  times  indicate  Moho  depths, 
upper  mantle  velocities,  and  velocity  gradients  with  values  in  the  south  and  the  north 
similar  to  those  found  for  event  Lie.  The  Pn  travel  times  are  consistent  with  the  a 
northward  shallowing  of  the  Moho  occurring  over  a  distance  of  roughly  150  km  between 
the  latitudes  31°  and  32.5°N,  similar  to  Lie. 

These  results  also  agree  with  the  data  from  events  Lli,  Llj,  and  Llk  with  epicenters 
located  near  event  Lla  and  are  shown  in  Figures  5(c),  (d)  and  (e).  Since  these  events 
have  similar  epicentral  locations,  they  provide  a  good  assessment  of  the  overall  data 
quality.  The  value  of  RMS  misfit  between  the  observed  travel  times  and  the  calculated 
travel  times  for  event  Lla  using  the  3D  regional  velocity  model  (Figures  6  and  7)  is  0.16 
s  and  similar  values  were  found  for  events  Lli,  Llj,  and  Llk. 

Event  E3b  is  located  to  the  northeast  of  the  Hi-CLIMB  array  and  Figure  5(f) 
shows  the  observed  travel  times  for  this  event.  The  calculated  Pg  travel  times  agree  with 
the  observed  travel  times  using  the  refined  crustal  velocity  model  determined  from  the 
four  events  located  nearby  the  array.  The  observed  Pn  travel  times  indicate  Moho  depths, 
upper  mantle  velocities,  and  velocity  gradients  with  values  in  the  south  and  north  parts  of 
the  model  similar  to  those  found  for  the  other  far  regional  events.  However,  the 
transition  between  Pn  travel  times  in  the  north  and  in  the  south,  shown  by  the  bar  in 
Figure  5(f),  indicates  that  the  Moho  shallows  northward  over  a  slightly  longer  distance  of 
around  200  km  and  slightly  farther  to  the  south  between  the  latitudes  of  30°  and  32°N. 
The  regional  velocity  model  determined  in  this  study  (Figures  6  and  7)  results  in  an  RMS 
misfit  between  the  calculated  and  observed  travel  times  for  this  event  of  0.26  s.  Some  of 
the  misfit  could  result  from  additional  variations  in  the  Pn  speeds  away  from  the  Hi- 
CLIMB  array,  and  also  in  the  vicinity  of  the  transition  in  Moho  depth  near  the  BNS. 

Results  and  Discussion 

Using  3D  ray  tracing  (Cerveny  et  al.,  1988;  the  CRT  algorithm)  to  forward  model 
the  Pg  and  Pn  travel  times,  we  have  determined  a  3D  P- wave  velocity  model  for  the  crust 
and  upper  mantle  in  the  central  and  western  Tibetan  Plateau.  Figure  6(a)  illustrates  the 
crustal  velocities  found  beneath  the  Hi-CLIMB  array  which  varies  laterally  from  north  to 
south  with  slightly  lower  velocities  in  the  crust  observed  in  the  areas  of  the  IYS  and  north 
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of  the  BNS.  Figure  7  highlights  the  velocity  structure  of  the  upper  mantle  and  shows  that 
upper  mantle  velocities  of  around  8.3  km/s  and  crustal  thicknesses  of  close  to  73  km  are 
found  beneath  the  southern  portion  of  the  Tibetan  Plateau.  To  the  north  of  the  BNS  in  the 
Qiangtang  terrane,  upper  mantle  velocities  are  found  to  be  lower  at  around  7.8  km/s,  and 
the  crust  in  this  area  is  about  10  km  thinner  than  that  observed  to  the  south. 

The  separation  between  the  Lhasa  terrane  to  the  south  and  the  Qiangtang  terrane  to 
the  north  can  also  be  observed  in  receiver  function  images  from  Hi-CLIMB  (Nabelek  et 
al.,  2009;  Nowack  et  al.,  2010)  as  a  disruption  in  the  Moho  conversion  in  the  receiver 
function  images  to  the  south  of  the  BNS.  The  location  of  this  disruption  coincides  with 
the  transition  in  the  Moho  depth  and  upper  mantle  velocities  found  in  this  study.  Figure 
6(b)  shows  a  comparison  of  the  Moho  depth  results  in  this  study  (gray  line)  from  wide- 
angle  travel  time  modeling  using  regional  events  with  the  teleseismic  imaging  results  of 
Nowack  et  al.  (2010),  where  the  solid  lines  are  the  interpreted  locations  of  the  subsurface 
Moho  and  other  structures.  Tseng  et  al.  (2009)  also  determined  Moho  depths  from  virtual 
wide-angle  seismic  profiles  using  teleseismic  data  and  their  wide-angle  imaging  results 
are  consistent  with  the  results  of  this  study  using  regional  wide-angle  events.  While  the 
results  from  geologic  mapping  and  balanced  cross-sections  given  in  Kapp  et  al.  (2005) 
and  Ding  et  al.  (2007)  suggest  a  northward  dipping  suture  near  the  surface  between  the 
Lhasa  and  Qiangtang  terranes,  the  results  of  this  and  other  geophysical  studies  from  Hi- 
CLIMB  (Tseng  et  al.,  2009;  Nowack  et  al.,  2010)  find  a  zone  of  changing  Moho  depth 
which  suggests  a  geophysical  suture  at  depth  mainly  located  to  the  south  of  the  surface 
expression  of  the  BNS  which  could  imply  a  southward  dipping  suture  between  the  Lhasa 
and  the  Qiangtang  terranes. 

Regional  tomography  studies  have  found  lower  than  average  Pg  and  Sg  velocities 
in  the  Tibetan  crust  of  around  6.1  and  3.4  km/s,  respectively  (Steck  et  al.,  2009; 

Monsalve  et  al.,  2008;  Sun  and  Toksoz,  2006).  Prior  studies  using  both  Hi-CLIMB  and 
INDEPTH  III  data  have  also  observed  low  velocity  zones  in  the  crust  near  the  areas  of 
the  IYS  and  BNS  (Hung  et  al.,  2010;  Meissner  et  al.,  2004),  which  this  is  similar  to  the 
results  found  from  this  study  for  the  crust.  Also,  the  results  for  upper  mantle  Pn 
velocities  found  in  this  study  are  consistent  with  earlier  studies  of  Pn  in  Tibet,  including 
those  of  Hearn  et  al.  (2004),  Sun  and  Toksoz  (2006),  Liang  and  Song  (2006),  and  Phillips 
etal.  (2007). 

The  inferred  upper  mantle  velocities  and  gradients  from  this  study  are  consistent 
with  those  of  Phillips  et  al.  (2007)  in  Tibet  with  upper  mantle  velocity  gradients  in  this 
region  of  Tibet  between  0.003  s'1  to  0.004  s'1.  Myers  et  al.  (2010)  found  similar  upper 
mantle  velocity  gradients  in  southern  Tibet,  but  lower  gradients  in  the  Qiangtang  terrane 
in  northern  Tibet.  However,  their  study  was  much  larger  in  scale  incorporating  Eurasia 
and  North  Africa.  For  this  study,  the  upper  mantle  velocities  are  well  modeled  by  the 
regional  travel  time  data,  but  the  modeled  upper  mantle  velocity  gradients  are  less 
constrained  (Griffin,  2010).  For  a  better  resolution  of  upper  mantle  velocity  gradients, 
longer  offsets  than  those  used  here  would  be  needed,  since  the  offsets  in  this  study  were 
restricted  to  include  only  events  with  well  recorded  cross-over  distances  between  the  Pg 
and  Pn  branches.  Alternatively,  other  seismic  attributes,  such  as  Pn  amplitudes,  could  be 
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used  to  provide  better  constraints  on  the  upper  mantle  velocity  gradients  and  is  the 
subject  of  current  research. 

Although  the  event  coverage  in  this  study  is  limited,  In  addition  to  the  north-south 
variation,  the  inferred  velocity  model  also  shows  an  east-west  lateral  variation.  The 
location  and  width  of  the  transition  in  Moho  depth  is  largely  constrained  by  the  travel 
times  for  ray  paths  from  the  regional  event  Lie  sampling  mostly  to  the  west  of  the  Hi- 
CLIMB  array  and  the  travel  times  of  the  ray  paths  of  event  E3b  sampling  predominantly 
to  the  east  of  the  array.  To  the  west,  the  Moho  is  found  to  shallow  to  the  north  between 
the  latitudes  of  31°  to  32.5°N,  whereas  to  the  east,  the  northward  shallowing  of  the  Moho 
is  found  to  occur  farther  to  the  south  between  30  to  32°N  and  generally  follows  the 
surface  trend  of  the  Bangong-Nujiang  suture  (BNS). 

Conclusions 

This  study  uses  broadband  seismic  data  recorded  by  the  Hi-CLIMB  array  in 
central  and  western  Tibet  from  earthquake  sources  located  both  along  and  at  more 
regional  distances  from  the  Hi-CLIMB  array.  These  events  provide  travel  time  data 
along  the  profile  of  the  array  but  also  afford  some  azimuthal  coverage  of  the  region. 
Observed  f-wave  travel  times  from  selected  regional  earthquakes  were  modeled  using 
3D  ray-tracing  in  order  to  determine  the  velocity  structure  of  the  crust  and  upper  mantle 
structure  in  Tibet  near  the  Hi-CLIMB  array.  Although  only  a  limited  number  of  events 
were  used  for  the  travel  time  modeling,  the  results  are  consistent  with  earlier  results  from 
teleseismic  imaging  beneath  the  Hi-CLIMB  array. 

A  laterally  varying  refined  crustal  velocity  model  was  found  using  Pg  travel  times 
from  the  events  located  along  and  nearby  the  Hi-CLIMB  array,  and  the  inferred  velocities 
are  similar  to  those  given  by  Meissner  et  al.  (2004)  for  INDEPTH  III  in  the  upper  crust 
but  with  slower  velocities  in  the  mid  to  lower  crust  beneath  the  Hi-CLIMB  array.  Crustal 
velocities  in  the  top  30  km  of  the  crust  were  also  found  to  be  up  to  0. 1 5  km/s  slower  in 
the  area  of  the  Indus-Y arlung  suture  south  of  latitude  3 1  °N  and  in  the  area  of  the 
Bangong-Nujiang  suture  north  of  latitude  33.5°N. 

The  Pn  travel  times  observed  for  the  more  regional  events  to  the  Hi-CLIMB  array 
were  used  to  constrain  the  Moho  structure  and  upper  mantle  velocities  in  the  region.  It 
was  found  that  the  crust  is  as  thick  as  73  km  beneath  the  southern  part  of  the  plateau  in 
the  Lhasa  terrane  and  about  10  km  thinner  beneath  the  northern  area  of  the  plateau  in  the 
Qiangtang  terrane.  This  north-south  variation  is  also  present  in  the  inferred  variation  in 
upper  mantle  velocities  which  were  found  to  be  about  8.3  km/s  beneath  the  Lhasa  terrane 
in  the  southern  plateau,  but  as  low  as  7.8  km/s  to  the  north  of  the  BNS  beneath  the 
Qiangtang  terrane.  However  the  travel  time  data  from  events  Lie  to  the  west  of  the  Hi- 
CLIMB  array,  and  E3b  to  the  east  of  the  array,  imply  variations  in  the  east-west  direction 
as  well,  generally  following  the  surface  trend  of  the  BNS.  The  transition  from  faster 
upper  mantle  velocities  and  a  thicker  crust  in  the  south  to  slower  upper  mantle  velocities 
and  a  thinner  crust  in  the  north  occurs  over  a  distance  of  150  km  between  the  latitudes  of 
31°  and  32.5°  N  in  the  western  portion  of  the  velocity  model  and  occurs  over  a  distance 
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of  around  200  km  between  the  latitudes  of  30°  and  32°N  in  the  eastern  portion  of  the 
model. 

Data  and  Resources 

Seismic  data  used  in  this  study  were  collected  as  part  of  the  Hi-CLIMB  Project 
using  instruments  from  a  number  of  different  sources,  including  many  from  the 
Incorporated  Research  Institutions  in  Seismology-Program  for  the  Array  Seismic  Studies 
of  the  Continental  Lithosphere  (IRIS-PASSCAL)  Center.  The  project  was  organized  and 
coordinated  jointly  by  the  University  of  Illinois  (principal  investigator,  Wang-Ping  Chen) 
and  Oregon  State  University  (principal  investigator,  John  Nabelek).  Seismic  data  can  be 
obtained  directly  from  the  IRIS  Data  Management  Center  at  www.iris.edu. 

Ray  tracing  was  performed  using  the  3D  ray  tracing  code  by  Cerveny  et  al.  (1988; 
the  CRT  algorithm),  and  event  relocation  was  performed  using  the  location  code 
HYPOINVERSE-2000  (Klein,  2002). 

The  USGS/NEIC  (PDE)  earthquake  locations  are  available  at 
http : //earthquake .usgs . gov/earthquakes/ eqarchives/ epic/ database .php  (last  accessed, 
November  2010).  The  EHB  catalog  of  earthquakes  (Engdahl  et  al.,  1998)  is  available  at 
http://www.isc.ac.uk/search/bulletin/ehb.html  (last  accessed,  October  2010). 
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Figure  Captions 


Figure  1.  The  topography  of  the  Tibetan  Plateau  shown  with  the  stations  of  the  Hi- 
CLIMB  array  in  Tibet  (triangles)  and  the  epicenters  of  the  selected  regional  earthquake 
sources  (circles).  Black  circles  indicate  the  USGS/NEIC  PDE  and  EHB  reported 
hypocenters,  and  the  white  dots  indicate  the  epicenters  of  the  four  events  nearest  to  the 
array  which  have  been  relocated  in  this  study.  The  locations  of  the  Bangong-Nujiang 
suture  (BNS)  and  the  Indus- Yarlung  suture  (IYS)  are  shown  by  the  dotted  lines,  and  the 
solid  line  gives  the  reference  line  along  the  Hi-CLIMB  array  used  for  plotting  the 
locations  of  the  events  L2a,  L2d,  L3a,  and  L3b  as  used  in  Figure  4. 

Figure  2.  Examples  of  vertical  component  record  section  recorded  by  the  Hi-CLIMB 
array  in  central  Tibet,  (a)  Shows  the  record  section  for  earthquake  L3a,  and  (b)  shows  the 
record  section  for  earthquake  Lie.  On  both  record  sections,  the  first  arriving  Pg  or  Pn 
waves  are  indicated.  The  waveforms  were  filtered  with  a  six  pole  Butterworth  bandpass 
filter  between  .5  and  5.0  Hz.  The  data  are  shown  in  reduced  time  by  8.0  km/s,  and  each 
trace  is  amplitude  normalized.  For  each  trace  on  the  record  sections,  the  first  arrivals  of 
P-waves  are  marked  by  black  dots. 

Figure  3.  A  rendering  of  the  P- wave  ray-paths  calculated  by  the  3D  ray-tracing  code 
CRT  (black  lines)  for  the  event  Lie  (white  sphere)  as  viewed  from  the  west.  The  gray 
surface  indicates  the  structure  of  the  Moho  used  within  the  input  velocity  model  for  the 
ray  tracing. 

Figure  4.  (a)  The  initial  ID  crustal  velocity  model  (solid  line)  found  for  the  region  of  Hi- 
CLIMB  with  upper  crustal  velocities  similar  to  those  found  by  Meissner  et  al.  (2004)  but 
with  slower  velocities  in  the  mid  to  lower  crust,  (b)  The  observed  travel  times  of  the  first 
arriving  P-waves  for  the  relocated  events  L2a  (diamonds),  L2d  (squares),  L3a  (triangles), 
and  L3b  (circles)  compared  with  the  calculated  travel  times  (lines  and  crosses)  using  the 
initial  ID  crustal  velocity  model.  All  travel  times  are  reduced  by  8.0  km/s  and  plotted 
with  respect  to  the  minimum  recorded  travel  time.  The  origin  of  each  event  is  placed 
respective  to  its  projected  position  along  the  solid  line  in  Figure  1.  (c)  The  same 
observations  as  in  (b)  but  with  travel  times  (lines  and  crosses)  predicted  from  using  the 
refined  crustal  velocity  model  and  laterally  varying  Moho  model  shown  in  Figure  6  and 
7. 

Figure  5.  The  observed  P-  wave  travel  times  (circles)  for  events  (a)  Lie,  (b)  LI  a,  (c)  Lli, 
(d)  Llj,  (e)  Llk,  and  (f)  E3b.  These  sources  are  located  at  regional  distances  to  the  Hi- 
CLIMB  array,  and  Pn  travel  times  for  these  events  were  used  to  constrain  the  P- wave 
velocity  structure  of  the  upper  mantle  for  the  region  of  the  Tibetan  Plateau.  All  travel 
times  are  reduced  by  8.0  km/s  and  plotted  with  respect  to  the  minimum  recorded  travel 
time.  The  lines  and  crosses  in  each  figure  show  the  calculated  travel  times  using  the  3D 
regional  velocity  model  found  in  this  study  and  shown  in  Figs.  6  and  7.  Horizontal  bars 
highlight  the  locations  of  the  transition  for  shoaling  of  the  Moho  and  northward  decrease 
in  Pn  speed. 
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Figure  6.  (a)  The  C- velocity  model  determined  in  this  study  showing  a  2D  transect  along 
the  strike  of  the  Hi-CLIMB  array  shown  in  Figure  1 .  Solid  lines  show  velocity  contours 
in  km/s.  Generally,  slower  crustal  velocities  were  found  beneath  the  areas  of  the  BNS 
and  IYS.  Also,  a  thicker  crust  and  faster  upper  mantle  velocities  were  found  in  the 
southern  area  of  the  plateau  beneath  the  Lhasa  terrane.  (b)  Comparison  of  the  crustal 
thickness  beneath  the  area  of  the  Hi-CLIMB  array  determined  in  this  study  (gray  line) 
with  the  interpretation  of  the  major  impedance  discontinuities  (black  lines)  adapted  from 
the  teleseismic  imaging  by  Nowack  et  al.  (2010). 

Figure  7.  A  rendering  of  the  3D  Moho  and  upper  mantle  velocity  structure  found  in  this 
study  for  the  region  of  the  Tibetan  Plateau.  The  upper  mantle  velocities  are  higher  and 
the  Moho  deeper  beneath  the  southern  portion  of  the  plateau  south  of  the  BNS,  and  upper 
mantle  velocities  are  lower  and  the  Moho  is  shallower  north  of  the  BNS.  However,  the 
area  of  higher  velocities  extends  farther  to  the  north  in  the  west  than  in  the  east  and 
generally  follows  the  trend  of  the  Bangong-Nujiang  suture  (BNS)  as  shown  by  the  dashed 
lines.  The  Indus-Yarlung  suture  (IYS)  is  also  shown. 
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FIGURE  4 
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Event  Lie  Compared  with  3D  Model 
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FIGURE  6 


70 


FIGURE  7 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


AFRL 

Air  Force  Research  Laboratory 

BNS 

Bangong-Nuijiang  Suture  in  Tibet 

CMT 

Centroid  moment  tensor  catalog 

CRT 

Complete  Ray  Tracing  computer  package 

MPI 

Message  Passing  Interface 

PDE 

Preliminary  Determination  of  Epicenters  catalog 

SEM 

Spectral  element  method 
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